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ABSTRACT 

This  thesis  introduces  a  new  procedure  for  the  enhancement  of  acoustic  images  of  the 
bottom  of  the  sea  produced  by  side-scan  sonars.  Specifically,  it  addresses  the  problem  of 
estimating  and  correcting  geometric  distortions  frequently  observed  in  such  images  as  a 
consequence  of  motion  instabilities  of  the  sonar  array.  This  procedure  estimates  the  ge¬ 
ometric  distortions  from  the  image  itself,  without  requiring  any  navigational  or  attitude 
measurements.  A  mathematical  model  for  the  distortions  is  derived  from  the  geometry 
of  the  problem,  and  is  applied  to  estimates  of  the  local  degree  of  geometric  distortion 
obtained  by  cross-correlating  segments  of  adjacent  lines  of  the  image.  The  model  param¬ 
eters  are  then  recursively  estimated  through  deterministic  least-squares  estimation.  An 
alternative  approach  based  on  adaptive  Kalman  filtering  is  also  proposed,  providing  a 
natural  framework  in  which  a  priori  information  about  the  array  dynamics  may  be  easily 
incorporated.  The  estimates  of  the  parameters  of  the  distortion  model  are  used  to  rectify 
the  image,  and  may  also  be  used  for  estimating  the  attitude  parameters  of  the  array.  A 
simulation  is  employed  to  evaluate  the  effectiveness  of  this  technique  and  examples  of  its 
application  to  high-resolution  side-scan  sonar  images  are  provided. 
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Chapter  1 


Introduction 


This  thesis  presents  a  new  procedure  for  the  enhancement  of  acoustic  images  made  by 
side-scan  sonars  in  the  sea  and  other  bodies  of  water.  Specifically,  a  solution  is  proposed 
to  the  problem  of  correcting  geometric  distortions  often  observed  in  those  images  as  a 
result  of  sonar  motion  instabilities.  The  techniques  introduced  here  are  unique  in  that 
they  require  no  navigational  data  or  attitude  measurements  for  the  sonar  array  to  correct 
the  geometric  distortions  in  the  image. 

The  remainder  of  this  chapter  presents  a  more  detailed  description  of  the  subject  of 
this  thesis.  An  introduction  to  side-scan  sonar  is  given  in  Section  1.1,  including  a  brief 
review  of  the  history  of  its  development,  its  principles  of  operation  and  applications. 
An  example  is  shown  of  a  side-scan  sonar  image  taken  from  our  data  set.  Section  1.2 
gives  a  description  of  the  various  types  of  distortions  frequently  observed  in  side-scan 
sonar  records,  with  emphasis  on  the  geometric  distortions  on  which  this  thesis  is  focused. 
The  current  state  of  research  on  digital  processing  of  side-scan  sonar  data  is  reviewed 
in  Section  1.3,  providing  a  framework  in  which  the  contributions  of  this  thesis  can  be 
situated.  Section  1.4  states  the  goals  of  our  work.  Finally,  Section  1.5  describes  the 
survey  during  which  the  data  was  collected. 

Chapter  2  presents  the  theoretical  development  of  a  mathematical  model  relating  ge¬ 
ometric  distortions  in  the  images  to  fluctuations  in  the  sonar  trajectory  and  heading. 
The  measurement  of  the  local  degree  of  geometric  distortion  in  small  areas  of  the  image 
from  the  cross-correlation  of  line  segments  is  discussed  in  Chapter  3.  Chapter  4  treats 
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the  estimation  of  distortion  parameters  using  the  model  derived  in  Chapter  2  and  the 
measurements  obtained  in  Chapter  3.  In  Chapter  5  the  parameter  estimates  are  used 
for  correcting  the  geometric  distortions  in  sonographs  and  examples  of  processed  images 
are  given.  A  simulation  is  used  for  evaluating  the  accuracy  of  the  geometric  correction. 
Chapter  6  presents  an  evaluation  of  the  results  obtained  in  this  thesis,  along  with  con¬ 
cluding  remarks  and  suggestions  for  future  research  on  this  problem.  Finally,  Appendix  A 
presents  a  summary  of  the  variables  defined  throughout  the  thesis,  and  Appendix  B  con¬ 
tains  the  derivation  of  some  of  the  results  from  the  geometrical  analysis  carried  out  in 
Chapter  2. 

1.1  An  Introduction  to  Side-Scan  Sonar 

The  development  of  diving  and  submersible  technology  has  made  the  bottom  of  the 
sea  a  stage  for  a  large  variety  of  enterprises.  Experiments  and  studies  in  marine  biol¬ 
ogy,  seabed  geology,  physical  oceanography,  and  other  scientific  disciplines  are  routinely 
carried  out  on  the  seafloor.  Archeologists  and  adventurers  alike  seek  and  explore  old 
ship  wrecks.  Communication  and  power  cables  are  laid  on  the  bottom  while  engineering 
structures  are  built  or  deployed  on  the  seabed  for  extracting  oil,  gas,  and  other  minerals, 
or  for  logistic  support  of  divers.  Many  of  these  activities  are  carried  out  not  only  in  the 
sea,  but  on  the  bottom  of  lakes  and  rivers  as  well. 

The  establishment  of  all  these  underwater  activities  has  generated  an  increasing  need 
for  underwater  remote  sensing.  This  is  often  provided  by  acoustic  bathymetric  devices, 
such  as  the  basic  echo  sounder  and  its  more  sophisticated  multi-beam  versions,  which  have 
been  extensively  used  to  produce  bathymetric  charts  of  the  oceans.  There  are,  however, 
many  underwater  enterprises  that  require  remote  sensing  at  a  much  finer  resolution  than 
can  be  achieved  with  bathymetric  devices.  In  some  cases,  visualization  of  the  bottom 
can  be  achieved  through  underwater  photography,  filming,  or  television.  However,  these 
means  are  severely  limited  by  the  rapid  absorption  of  light  by  the  water,  resulting  in 
maximum  ranges  of  at  most  a  few  tens  of  meters,  and  sometimes  as  little  as  less  than  a 
meter,  depending  on  the  condition  of  the  water. 
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Side-scan  sonar  bridges  the  gap  between  sight  and  sound  in  the  underwater  environ¬ 
ment  by  providing  visualization  of  the  bottom  through  acoustic  remote  sensing  at  ranges 
far  beyond  those  afforded  by  optical  means.  It  has  acquired  an  increasingly  important 
role  in  underwater  exploration  and  its  wide  range  of  applications  includes  mapping  the 
seabed  as  an  aid  to  the  production  of  nautical  charts  and  for  oceanographic  and  under¬ 
water  archeologic  research,  search  and  location  of  objects  on  the  bottom,  fishery  studies, 
support  for  submersible  operations,  and  off-shore  mining  and  engineering  surveys.  The 
history,  principles  of  operation,  and  applications  of  side-scan  sonars  are  reviewed  in  [58], 
[31],  [7],  [48],  [52],  and  [3]. 

1.1.1  History 

The  origins  of  side-scan  sonar  can  be  traced  back  to  the  submarine  research  conducted 
after  World  War  II  in  England,  where  sonar  operators  of  the  Allied  Submarine  Devices 
Investigation  Committee  (ASDIC)  observed  that  the  intensity  of  backscattering  of  high- 
frequency  sonar  pulses  from  the  bottom  appeared  to  present  a  consistent  correlation  with 
seabed  morphology  [24].  Kunze,  in  1957  [35],  and  Chesterman  et  al.,  in  1958  [9],  were 
the  first  to  explore  this  phenomenon  in  an  attempt  to  produce  a  representation  of  the 
seabed  topography.  By  1961  the  first  operational  side-scan  sonar  had  been  built  by  the 
National  Institute  of  Oceanography  (NIO)  in  England,  and  mounted  on  the  hull  of  its 
research  ship,  RRS  Discovery  II  [59].  The  first  major  seabed  survey  conducted  with  that 
sonar  was  reported  in  1961  [19]. 

The  development  of  side-scan  sonar  in  the  United  Kingdom  has  been  carried  out  by  the 
scientific  community,  especially  by  the  NIO,  now  known  as  the  Institute  of  Oceanographic 
Sciences  (IOS).  Their  most  notable  achievement  was  the  Geological  Long-Range  Inclined 
Asdic  (GLORIA)  Project1,  responsible  for  the  development  of  a  long-range  side-scan 
sonar  completed  in  1969  [47].  In  1977  the  IOS  built  a  second  version,  the  GLORIA  II, 
which  incorporated  a  series  of  improvements  suggested  by  the  experience  gained  in  the 
operation  of  the  original  GLORIA  sonar  [51]. 

JIn  its  early  days,  side-scan  sonar  was  sometimes  called  Asdic  sonar,  because  of  its  association  with 
the  Allied  Submarine  Devices  Investigation  Committee. 
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In  the  United  States,  the  development  of  side-scan  sonar  was  carried  out  not  only 
by  scientific  institutions  but  also  by  the  industry,  and  American  companies  have  held  a 
prominent  position  in  the  world  market  from  the  beginning.  Notable  examples  of  side- 
scan  sonars  developed  by  scientific  institutions  in  the  United  States  are  the  SeaMARC 
(Sea  Mapping  and  Remote  Characterization)  sonar  and  its  successor,  the  SeaMARC  II, 
jointly  developed  by  the  University  of  Hawaii  and  International  Submarine  Technology, 
Ltd.  [6] 

1.1.2  Principles  of  Operation 

Figure  1.1  depicts  a  typical  operational  setting  for  a  conventional  side-scan  sonar. 
The  main  components  of  the  system  are  the  so-called  towfish ,  a  cylindrical  towed  body 
on  whose  sides  are  mounted  two  linear  arrays  of  hydro-acoustic  transducers,  and  the 
recorder,  a  piece  of  equipment  containing  the  circuitry  that  generates  the  sonar  pulses 
and  detects  the  returned  signal,  producing  a  graphic  record  of  its  strength  either  on 
paper  or  on  a  cathode-ray  tube.  The  arrays  of  transducers  produce  two  sound  beams 
with  very  narrow  main  lobes  (typically  only  0.2  to  2.0  degrees  wide)  on  each  side  of  the 
towfish,  pointing  laterally  and  down  as  depicted  in  the  figure.  Each  sonar  pulse  ensonifies 
a  narrow'  strip  of  the  seabed  on  each  side,  to  a  distance  of  up  to  30  kilometers,  depending 
on  the  operating  frequency  of  the  sonar.  As  the  towfish  is  towed  through  the  water,  the 
lateral  sound  beams  progressively  scan  a  swath  of  the  seabed  (thus  the  name  side-scan 
sonar),  while  the  recorder  produces  a  line-by-line  record  of  the  backscattered  signal  [25]. 

The  visual  record  produced  by  a  side-scan  sonar  is  called  i  sonograph.  Figure  1.2 
shows  an  example  of  a  sonograph  covering  an  area  of  100  m  x  200  m  of  a  fairly  flat 
stretch  of  seabed  consisting  of  a  field  of  boulders  and  cables  laid  on  the  bottom.  The 
towfish  may  be  pictured  moving  along  the  center  of  the  image  from  bottom  to  top,  so 
that  the  right-hand  and  left-hand  halves  of  the  image  correspond  to  the  starboard  and 
port  sides,  respectively.  The  two  arrays  of  transducers  fire  a  sequence  of  sound  pulses  at 
regular  intervals  as  the  towfish  is  towed  above  the  seabed.  After  each  pulse  is  emitted, 
the  transducers  receive  the  signals  returned  from  the  bottom  on  the  starboard  and  port 
sides,  before  the  next  pulse  is  fired.  In  the  case  of  a  digitized  sonograph,  the  returned 
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Figure  1.1:  Typical  operational  setting  of  a  side-scan  sonar.  Lateral  sound  beams  scan 
the  seabed  as  the  sonar  towfish  is  towed  above  the  bottom  by  a  deploying  vessel. 

signals  are  sampled  at  a  constant  rate  to  form  one  line  of  the  image,  starting  from  its 
middle  point  and  moving  towards  the  sides.  The  image  intensity  is  proportional  to  the 
strength  of  the  returned  signal,  and  the  horizontal  coordinate  measured  from  the  middle 
of  the  line  corresponds  to  the  time  elapsed  since  the  pulse  was  sent  out.  Equivalently, 
the  horizontal  coordinate  may  be  viewed  as  range  from  the  towfish,  assuming  the  sound 
pulse  travels  at  a  constant  speed.  The  vertical  coordinate  of  the  sonograph  corresponds 
to  the  distance  traveled  by  the  towfish.  The  dark  stripe  seen  in  the  middle  of  the  image 
corresponds  to  the  interval  of  time  between  the  emission  of  the  sound  pulse  and  the 
instant  it  first  reaches  the  seabed,  during  which  there  are  no  reflections  except  those 
from  fish  or  other  objects  in  the  water  column. 

As  a  visual  record  of  backscattered  sound,  a  sonograph  depicts  a  combination  of  the 
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Figure  1.2:  Example  of  a  side-scan  sonar  image,  after  processing  for  correction  of  intensity 
distortions.  Geometric  distortions  are  responsible  for  the  jagged  appearance  of  the  cable 
seen  in  the  lower  right-hand  corner  and  for  multiple  images  of  several  rocks  in  that  area. 


seabed  topography  and  its  geological  constitution.  This  is  a  consequence  of  the  depen¬ 
dence  of  backscattering  strength  both  on  the  angle  of  incidence  at  which  the  sound  waves 
reach  the  target  (called  the  grazing  angle)  and  on  the  target’s  ability  to  reflect  them  back, 
that  is,  its  reflectivity  [11].  As  a  result,  sonographs  resemble  aerial  photographs,  with 
the  difference  that  the  coordinates  of  a  point  on  a  sonograph  correspond  to  horizontal 
distances  in  one  direction  (along-track)  and  to  range  to  the  towfish  in  the  other  direction 
(cross- track),  whereas  in  aerial  photographs  the  coordinates  of  a  point  in  both  direc¬ 
tions  correspond  to  horizontal  distances  on  the  ground.  A  closer  analogy  exists  between 
sonographs  and  side-looking  radar  images,  which  are  also  formed  line-by-line  and  whose 
coordinates  are  also  distance  and  range  [50]. 

1.1.3  Types  and  Applications  of  Side-Scan  Sonar 

The  most  relevant  parameter  in  the  characterization  of  a  side-scan  sonar  is  its  op¬ 
erating  frequency,  which  fixes  its  maximum  range  and,  together  with  the  array  length 
and  pulse  width,  determines  its  resolution.  Sound  absorption  in  water  increases  with 
frequency,  resulting  in  shorter  ranges  [11].  However,  the  beam  width  and  pulse  length 
can  be  made  smaller  at  higher  frequencies,  providing  increased  resolution  [60].  Thus, 
side-scan  sonars  with  low  operating  frequencies  (3.5  to  40  kHz)  have  very  long  ranges  (1 
to  30  km  on  each  side)  but  knver  resolution,  while  high-frequency  units  (100  to  500  kHz) 
have  resolutions  as  fine  as  a  few  centimeters  at  the  expense  of  shorter  maximum  ranges 
(100  m  to  1  km)  [24].  The  former  class  of  sonars  is  mostly  employed  for  marine  geolog¬ 
ical  and  mapping  surveys,  where  it  is  necessary  to  cover  very  large  areas  of  the  seabed 
and  only  the  large-scale  topography  is  of  interest.  The  latter  class  is  better  suited  to 
applications  such  as  the  visualization  of  objects  on  the  bottom,  which  require  increased 
resolution  over  smaller  areas. 

A  notable  example  of  a  long-range  side-scan  sonar  is  the  GLORIA  system,  which, 
with  its  long  total  coverage  of  up  to  60  km  (at  an  operating  frequency  of  6.5  kHz)  has 
produced  a  wealth  of  data  over  15  years  of  operation  by  the  IOS  and  on  charters  to 
foreign  organizations  [36].  The  United  States  Geological  Survey  has  employed  GLORIA 
for  mapping  extensive  portions  of  the  US  Exclusive  Economic  Zone  (EEZ).  An  example 
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of  a  high-frequency  unit  is  the  500  kHz  Klein  Hydroscan  sonar,  which,  with  its  very  fine 
resolution  of  less  than  5  cm,  has  produced  some  striking  images  of  shipwrecks  and  other 
objects  on  the  bottom  (21]. 

A  number  of  variations  have  evolved  from  the  conventional  side-scan  sonar.  Perhaps 
the  most  significant  of  these  is  the  bathymetric  side-scan  sonar,  which  uses  the  phase 
of  the  backscattered  signal  at  each  moment  for  calculating  the  seafloor  elevation  as  a 
function  of  range  for  each  line  of  the  sonograph  (15]  [16]  [17].  An  example  of  a  bathymetric 
side-scan  sonar  is  the  SeaMARC  II  [6].  Another  variation  is  the  interferometric  side-scan 
sonar  [33].  The  application  of  synthetic-aperture  techniques  to  side-scan  sonar  has  also 
been  proposed,  but  success  in  this  area  has  been  limited  [49]  [27]. 

1.2  Sonograph  Distortions 

The  distortions  that  affect  sonographs  are  of  two  kinds:  intensity  distortions,  which 
are  deviations  from  the  ideal  linear  relation  between  image  intensity  and  backscattering 
strength  of  the  materials  on  the  bottom,2  and  geometric  distortions,  which  correspond 
to  discrepancies  between  the  relative  location  of  features  on  the  image  and  their  true 
location  on  the  seabed.  Various  specific  forms  of  distortions  can  be  identified  in  each  of 
these  two  categories  [56]  [23]  [8].  The  different  types  of  distortion  will  be  discussed  in 
the  remainder  of  this  section. 

1.2.1  Intensity  Distortions 

Several  types  of  intensity  distortions  can  be  observed  in  sonographs.  One  of  them  is 
caused  by  power  drop-off  in  the  returned  signal,  resulting  from  three  factors:  The  first 
two  are  spherical  spreading  loss  and  sound  absorption  by  the  water,  both  of  which  cause 
the  sound  wave  to  be  attenuated  as  it  travels  away  from  the  towfish  and  back.  The  third 
factor  is  the  decrease  in  backscattering  due  to  diminishing  grazing  angles  at  increasingly 
longer  ranges.  Most  side-scan  sonars  attempt  to  correct  power  drop-off  by  applying  a 

3These  are  sometimes  called  niiomeiric  distortions,  a  term  apparently  borrowed  from  the  radar 
remote  sensing  nomenclature. 
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time-varying  gain  to  the  returned  signal,  which  the  user  can  adjust  to  obtain  an  image 
intensity  as  uniform  as  possible.  In  the  process,  one  may  eliminate  some  legitimate  large- 
scale  intensity  variations  caused  by  a  sloping  bottom  or  by  differences  in  the  geological 
constitution  of  the  seabed. 

A  source  of  intensity  distortions  inherent  to  the  sonar  itself  is  the  effect  of  the  side 
lobes  of  the  beam  pattern.  There  are  both  vertical  sidelobes,  i.  e.,  those  on  the  plane 
perpendicular  to  the  axis  of  the  towfish,  and  horizontal  sidelobes,  i.  e.,  those  outside  that 
plane.  The  vertical  sidelobes  beneath  the  main  lobe  may  be  responsible  for  variations  in 
intensity  in  the  first  portion  of  the  returned  signal,  while  those  above  the  main  lobe  may 
cause  spurious  surface  reflections.  In  the  cross-track  direction,  each  line  of  the  sonograph 
is  produced  by  the  convolution  of  the  acoustic  beam  function,  defined  by  the  main  lobe 
and  horizontal  sidelobes,  with  the  backscattering  function  of  the  seabed.  As  a  result, 
each  line  is  not  the  perfect  sample  of  the  bottom  that  one  would  desire.  This  effect  is 
usually  more  noticeable  at  far  ranges,  where  the  beam  spreading  may  cause  the  features 
to  appear  smeared  in  the  along-track  direction.  Pitching  and  yawing  of  the  towfish  also 
result  in  variations  in  image  intensity  by  causing  misalignment  of  the  transmitting  and 
receiving  beams,  with  a  consequent  drop  in  the  power  level  of  the  returned  signal,  which 
translates  into  striping  of  the  sonograph  in  the  cross-track  direction. 

Other  types  of  intensity  distortions  are  ultrasonic  interference  generated  by  passing 
ships,  electrical  interference  from  other  instruments  on  the  ship,  and  various  kinds  of 
artifacts  caused  by  shoals  of  fish,  dense  particle  suspension  in  the  water,  or  other  distur¬ 
bances  of  the  medium. 

1.2.2  Geometric  Distortions 

One  type  of  geometric  distortion  is  inherent  to  side-scan  sonar:  the  so-called  slant- 
range  distortion ,  a  consequence  of  the  cross-track  coordinate  of  sonographs  being  range 
to  the  towfish  rather  than  horizontal  distance  on  the  bottom.  This  distortion  amounts 
to  a  range-dependent  coordinate  transformation  that  compresses  the  scale  towards  the 
shorter  ranges  and  causes  the  water  column  to  appear  on  the  sonograph  as  depicted  in 
Fig.  1.3.  If  the  bottom  is  assumed  to  be  flat,  correcting  the  slant-range  distortion  is  a 
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Figure  1.3:  The  slant-range  distortion  is  caused  by  the  fact  that  horizontal  distances  in 
the  sonograph  correspond  to  ranges  to  the  towfish,  rather  than  distances  on  the  seabed. 

fairly  simple  operation  that  only  requires  knowledge  of  the  height  of  the  towfish  above 
the  bottom.  That  height  can  be  measured  from  the  sonograph  itself  as  the  width  of  the 
water  column  at  each  line,  and  then  be  used  to  convert  ranges  to  horizontal  distances  on 
the  bottom.  This  procedure  will  be  described  in  greater  detail  in  Section  6.1. 

Another  source  of  geometric  distortions  is  variations  of  the  speed  of  sound  in  the  water, 
caused  by  differences  in  water  temperature,  pressure,  or  salinity  [11].  Gradients  in  the 
velocity  of  sound  in  the  water  also  cause  the  sound  waves  to  be  refracted,  a  phenomenon 
known  as  ray  bending  [18].  Both  these  phenomena  preclude  an  exact  correspondence 
between  sound-wave  travel  time  and  actual  range  to  the  towfish,  with  the  result  that 
cross-track  distances  on  the  sonograph  no  longer  represent  actual  ranges.  Under  extreme 
conditions,  ray  bending  may  cause  the  beam  to  miss  some  target  areas  altogether,  leaving 
shadow  zones  in  the  sonograph.  These  kinds  of  distortion  are  usually  observed  only  at 
longer  ranges  and  are  uncommon  in  high-frequency  sonographs. 

Still  another  source  of  geometric  distortion  in  sonographs  is  variations  in  the  trajec- 
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tory,  speed,  or  orientation  of  the  towfish.  Ideally,  the  towfish  would  be  towed  above  the 
bottom  at  a  constant  speed,  on  a  straight  path,  and  with  the  heading  always  aligned  with 
the  trajectory.  In  practice,  however,  the  towfish  is  often  subject  to  motion  instabilities. 
Wind  and  sea  currents,  for  instance,  often  prevent  the  deploying  vessel  from  maintain¬ 
ing  a  constant  speed  and  heading.  As  a  result,  the  sonograph  will  display  variations  in 
aspect-ratio  and  other  large-scale  geometric  distortions.  These  can  be  corrected  by  using 
navigational  measurements  to  determine  the  true  location  of  the  ship  and  towfish  and 
then  resampling  the  sonograph  so  as  to  produce  an  isometric  image.  Towfish  motion 
instabilities  may  also  result  from  the  action  of  underwater  currents  on  the  towfish  or 
from  the  ship  sway  being  communicated  to  the  towfish  through  the  tow  cable. 

Motion  instabilities  may  be  divided  into  two  types:  translational  and  rotational. 
Translational  instabilities  correspond  to  lateral  and  vertical  displacements  of  the  towfish 
from  the  desired  straight  path,  and  to  variations  in  speed.  Rotational  instabilities  are 
illustrated  in  Fig.  1.4.  Pitching  and  yawing  produce  geometric  distortions  by  causing  the 
beams  to  scan  ahead  or  back,  either  simultaneously  on  both  sides,  in  the  case  of  pitching, 
or  alternating  between  the  port  and  starboard  sides,  in  the  case  of  yawing.  It  is  common 
for  pitching  and  yawing  to  be  severe  enough  to  cause  the  sonar  beam  to  scan  backwards 
over  a  previously  covered  area.  Objects  in  a  back-scanned  area  appear  in  triplicate  on  the 
sonograph:  one  image  corresponding  to  the  first  time  they  were  scanned,  followed  by  their 
mirror  image  reversed  in  the  along-track  direction  resulting  from  the  back-scanning,  and, 
finally,  a  third  image  produced  when  the  beam  starts  scanning  forward  again.3  Rolling 
does  not  produce  geometric  distortions,  but  it  causes  variations  in  image  intensity  in  the 
area  immediately  below  the  towfish  due  to  the  rotation  of  vertical  sidelobes.  It  can  also 
result  in  the  near-field  portion  of  one  side  of  the  sonograph  being  folded  back  on  itself 
by  allowing  the  lower  part  of  the  corresponding  beam  to  cross  over  into  the  opposite  side 
under  the  towfish. 


3Back-scanning  does  not  occur  in  the  cross-track  direction,  since  it  is  impossible  for  the  propagating 
wavefront  to  reverse  its  course  as  a  result  of  towfish  instabilities. 
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Figure  1.4:  Rotational  instabilities  of  the  towfish:  pitching  (top),  yawing  (middle),  and 
rolling  (bottom). 

1.3  Digital  Processing  of  Sonographs 

The  various  types  of  distortions  described  in  Section  1.2  can  sometimes  significantly 
limit  the  accuracy  of  sonographs  as  visual  representations  of  the  seabed.  It  was  not  long 
after  the  invention  of  side-scan  sonar  that  the  first  attempts  were  made  at  correcting 
sonograph  distortions.  For  about  a  decade,  the  only  types  of  distortions  that  were  cor¬ 
rected  were  those  caused  by  variations  in  ship  speed  and  by  the  slant-range  effect.  In  [10], 
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for  instance,  a  system  was  described  in  which  the  slant-range  distortion  w as  removed  by 
using  a  nonlinear  helix  in  the  recorder,  and  variations  in  ship  speed  were  compensated  for 
by  using  a  continuous  flow  camera  to  change  the  along-track  scale  of  the  sonograph.  The 
same  results  were  accomplished  by  fiber  optic  recorders  described  in  [53]  and  [61],  and  by 
a  system  reported  in  [4]  that  employed  a  photo-electrical  scanner  in  conjunction  with  a 
recorder  whose  stylus  was  driven  by  a  stepping  motor  in  a  non-linear  sweep.  Commercial 
systems  that  correct  for  slant  range  and  for  variations  in  ship  speed  were  later  developed, 
such  as  the  ones  described  in  [13]  and  [32]. 

Apparently,  the  first  system  reported  in  the  literature  that  departed  from  the  tradi¬ 
tional  method  of  recording  the  sonar  signal  on  paper  was  the  one  described  in  [28]  and 
[29].  In  that  system  the  sonar  signal  was  recorded  on  magnetic  tape  during  the  survey 
and  was  later  displayed  on  a  cathode-ray  tube,  affording  a  higher  dynamic  range  than 
the  paper  record  and  at  the  same  time  allowing  the  operator  to  correct  for  variations  in 
ship  speed  and  heading  by  manually  adjusting  the  vertical  time  base  and  the  slope  of  the 
scanning  lines. 

From  the  beginning,  recorders  have  incorporated  a  time- varying  gain  (TVG),  which 
is  applied  to  the  sonar  signal  in  an  attempt  to  compensate  for  the  power  drop-off  with 
range  caused  by  sound  attenuation,  spreading  loss,  and  decreasing  grazing  angles,  as 
explained  in  Section  1.2.1.  Beyond  that,  the  only  other  instance  of  analog  processing 
of  sonographs  appears  to  have  been  the  use  of  analog  filtering  for  selective  textural 
enhancement  reported  in  [34]. 

Correcting  for  these  simpler  geometric  and  intensity  distortions  was  all  that  could  be 
accomplished  without  the  recourse  of  digital  processing.  The  first  steps  in  that  direction 
were  taken  in  1976,  when  reports  appeared  in  the  literature  concerning  the  digitizing  of 
sonographs  [62]  and  their  processing  by  computers.  Paluzzi  et  al.  seem  to  have  been  the 
first  to  report  digital  processing  of  sonographs  [41],  Basically,  they  used  existing  software 
developed  at  the  Jet  Propulsion  Laboratory  of  the  California  Institute  of  Technology  for 
processing  images  from  unmanned  planetary  exploration  missions,  and  adapted  it  for 
sonograph  processing.  Use  of  navigational  data  for  correction  of  large-scale  geometric 
distortions  and  production  of  sonograph  mosaics  was  later  added  to  that  system,  and 
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it  was  employed  for  processing  GLORIA  sonographs  in  a  joint  project  with  the  British 
Institute  of  Oceanographic  Sciences  [42].  That  same  software  package  was  later  used 
and  extended  by  Luyendyk,  Hajic,  and  Simonett  at  the  University  of  California,  Santa 
Barbara  [39].  Similar  sonograph-processing  systems  have  since  been  developed  at  the 
National  Mapping  Division  of  the  U.  S.  Geological  Survey  as  reported  by  Chavez  in  [8] 
and  at  France’s  Institut  Frangais  de  Recherche  pour  l ’Exploitation  de  la  Mer  as  reported 
by  Augustin  in  [2].  The  systems  described  in  [22],  [20],  and  [46]  incorporate  techniques  for 
feature  extraction  and  image  segmentation  and  classification.  Other  instances  of  digital 
processing  of  sonographs  are  described  in  [12],  [45],  [49],  [57],  [55],  and  [30]. 

In  general,  one  might  say  that  considerable  improvements  are  still  possible  in  the 
field  of  digital  processing  for  sonograph  enhancement.  In  particular,  an  effective  tech¬ 
nique  for  correcting  the  geometric  distortions  caused  by  towfish  instabilities  without 
requiring  navigational  or  attitude  measurements  undoubtedly  constitutes  a  significant 
new  development. 


1.4  Objectives 

This  thesis  explores  the  problem  of  digitally  processing  sonographs  for  the  correction 
of  geometric  distortions.  Its  goal  is  to  introduce  techniques  for  estimating  and  correct¬ 
ing  these  distortions  using  solely  the  information  contained  in  the  image  itself,  without 
requiring  navigational  or  attitude  measurements. 

Geometric  distortions  in  the  digitized  images  are  interpreted  as  the  result  of  variations 
in  sampling  period  and,  from  the  geometry  of  the  problem,  a  model  is  derived  that  relates 
local  variations  in  the  sampling  period  to  parameters  that  describe  deviations  in  the 
trajectory  and  orientation  of  the  towfish.  A  linearized  version  of  the  model  is  then  used 
for  the  estimation  of  these  parameters. 

Measures  of  the  local  degree  of  geometric  distortion  of  the  image  are  derived  from 
the  cross-correlation  of  segments  of  adjacent  lines.  These  measurements  are  not  used 
directly  as  input  to  the  model,  but  rather  are  first  processed  to  yield  estimates  of  the 
sampling  intervals  on  the  seabed,  which  are  in  turn  used  for  the  estimation  of  distortion 
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parameters.  We  also  examine  practical  issues  concerning  the  calculation  of  the  sampling 
period  from  the  correlation  length,  such  as  the  detection  of  back-scanning,  and  the  effect 
of  the  sonar  beam  pattern. 

Estimation  of  the  distortion  parameters  from  the  correlation  length  is  accomplished 
through  deterministic  least-squares  estimation  on  a  line  by  line  basis.  An  alternative 
technique,  based  on  an  adaptive  Kalman  filtering  algorithm  known  as  the  extended  least- 
squares  method,  provides  a  framework  in  which  a  priori  knowledge  about  the  towfish 
dynamics  may  be  easily  incorporated  in  the  future.  Once  the  motion  parameters  are 
estimated,  the  geometric  distortions  are  corrected  by  resampling  the  image  appropriately. 


1.5  Description  of  the  Data  Set 

The  sonographs  used  in  this  thesis,  including  the  one  shown  in  Fig.  1.2,  were  collected 
during  a  survey  conducted  out  of  the  Woods  Hole  Oceanographic  Institution  (WIIOI) 
in  the  summer  of  1987.  Figure  1.5  shows  the  location  of  the  survey  area  off  the  island 
of  Martha’s  Vineyard,  Massachusetts.  That  area  was  chosen  because  of  the  existence  of 
utility  cables  on  the  bottom,  which  were  expected  to  facilitate  the  perception  of  geometric 
distortions  and  thus  provide  some  visual  indication  of  the  effectiveness  of  the  technique 
developed  for  correcting  them.  (The  technique  presented  in  this  thesis,  however,  does 
not  make  any  assumptions  about  or  rely  in  any  way  on  the  presence  of  such  objects  in 
the  image.) 

The  survey  was  conducted  aboard  WHOI’s  R.V.  Asterias  with  a  high-resolution  500- 
kHz  Klein  model  422S-001E  towfish  and  a  Klein  model  521  recorder  operating  at  a 
range  of  100  m,  the  maximum  afforded  by  that  sonar  unit.  The  demodulated  signal 
was  recorded  on  magnetic  tape  and  was  later  digitized  through  a  personal-computer 
system  equipped  with  a  Data  Translation  DT-2S51  “frame-grabber”  card.  The  size  of 
the  digitized  images  is  512  x  1024  pixels.  Assuming  a  sound  speed  of  1500  m/s,  the 
round-trip  time  of  the  sonar  pulse  corresponding  to  the  maximum  range  of  100  m  is 
133.33  ms,  which  determined  a  sampling  frequency  of  3.9  kHz  to  produce  512  pixels 
per  line  per  channel.  The  resulting  sampling  distance  in  the  across-track  direction  is 
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Figure  1.5:  Location  of  the  survey  conducted  in  the  summer  of  1987,  during  which  the 
sonographs  presented  in  this  thesis  were  collected. 

approximately  20  cm.  The  sonar  firing  rate  was  set  to  7.5  pulses/sec,  the  maximum 
possible  value  given  the  total  round-trip  travel  time  of  133.33  ms.  The  vessel  speed  was 
maintained  at  approximately  3  knots,  in  order  to  attain  the  same  sampling  distance  of 
approximately  20  cm  in  the  along-track  direction  as  in  the  cross-track  direction,  at  the 
selected  firing  rate. 

Though  the  sea  condition  during  the  survey  could  be  described  as  only  moderately 
rough,  the  small  size  of  the  vessel  (46  ft,  20  long  tons  loaded  displacement)  resulted  in 
enough  swaying  to  generate  a  significant  degree  of  geometric  distortion  in  the  images.  To 
that  was  added  the  effect  of  fairly  strong  undercurrents  prevailing  in  Vineyard  Sound. 
From  the  resulting  set  of  digitized  sonographs  a  number  of  images  were  chosen  to  be  used 
in  the  development  and  evaluation  of  the  techniques  described  in  this  thesis. 
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Chapter  2 


Mathematical  Analysis  of  the 
Problem 

This  chapter  presents  a  mathematical  analysis  of  geometric  distortions  in  sonographs, 
leading  to  the  development  of  a  linear  model  that  will  be  used  in  Chapter  4  for  estimating 
the  location  of  the  sampling  points  on  the  seabed.  The  analysis  presented  here  takes  into 
account  the  misalignment  of  the  transmitting  and  receiving  beams  of  the  sonar,  an  effect 
that  is  not  incorporated  in  previous  studies  found  in  the  side-scan  sonar  literature,  such 
as  [37],  [14],  and  [44]. 

2.1  Geometric  Distortions  in  Sonographs 

A  digitized  sonograph  can  be  viewed  as  a  mapping  of  the  backscattering  strength 
of  materials  on  the  seafloor  onto  a  monochrome  digital  image.  In  order  to  derive  a 
mathematical  description  of  this  transformation  it  is  necessary,  first  of  all,  to  intro¬ 
duce  a  three-dimensional  system  of  rectangular  coordinates  (z,y,z)  on  the  seabed  and  a 
two-dimensional  rectangular  coordinate  system  [m,n]  on  the  sonograph,1  as  depicted  in 
Fig.  2.1.  For  our  purposes,  the  backscattering  strength  of  the  bottom  can  be  completely 
described  by  a  function  b(x,y,z)  defined  on  the  seabed,  and  representing  the  backscat¬ 
tering  strength  at  point  (x,y,z).  As  discussed  in  Chapter  1,  the  backscattering  strength 

{The  square  brackets  are  used  to  indicate  that  m  and  n  are  discrete  rather  than  continuous  variables. 
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Figure  2.1:  Coordinate  systems  defined  on  the  seabed  and  on  the  sonograph. 

is  also  a  function  of  the  grazing  angle  at  which  the  towfish  detected  the  portion  of  the 
signal  reflected  by  point  ( x ,  y,  z)  and  of  the  operating  frequency  of  the  sonar,  although 
these  dependences  are  not  explicitly  indicated.  The  sonograph,  on  the  other  hand,  can 
be  described  by  a  two-dimensional  sequence  s[m,  n]  corresponding  to  the  image  intensity 
at  coordinates  [m,  n].  Here  m  denotes  the  image  coordinate  in  the  cross-track  direction, 
which  corresponds  to  lateral  ranges  to  the  towfish,  and  n  denotes  the  image  coordinate 
in  the  along-track  direction,  corresponding  to  the  time  elapsed  since  the  beginning  of 
the  acquisition  of  the  sonograph.  In  all  sonographs  shown  in  this  thesis,  the  direction  of 
scanning  is  from  bottom  to  top,  so  that  n  is  pointing  up,  while  m  is  positive  on  the  right 
or  starboard  side  and  negative  on  the  left  or  port  side  of  the  image. 

The  production  of  a  sonograph  is  then  denoted  by 

b[x,y,z)  *->  s[m,n]. 

This  operation  actually  comprises  an  intensity  transformation  b  *-»  3  and  a  sampling 
(z,y,  2)  *-*  [m,n].  The  intensity  transformation  is  not  a  simple  function  of  the  backscat- 
tering  coefficient,  being  affected  by  a  series  of  factors  such  as  sound  attenuation  and 
absorption  by  the  water,  the  angle  of  incidence  of  the  wavefront  on  the  bottom,  acoustic 
and  electrical  noise,  and  the  sonar  beam  pattern  and  time- varying  gain  [52],  making  it 
very  difficult  to  describe  mathematically.  The  sampling  operation,  on  the  other  hand,  is 
considerably  easier  to  describe  since  the  only  factors  that  it  involves  in  addition  to  the 
point  coordinates  are  the  position  and  orientation  of  the  towfish.  To  calculate  the  true 
location  of  points  of  the  sonograph  on  the  seabed,  it  is  necessary  to  determine  an  expres¬ 
sion  not  for  the  direct  sampling  operation,  but  for  the  inverse  mapping  [m,  n]  •-»  (x,  y,  z). 
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Before  proceeding,  we  make  a  simplifying  assumption.  Standard  side-scan  sonars  produce 
no  information  on  the  bearing  of  the  returned  signal,  and  it  is  impossible  to  determine 
all  three  of  x,  y,  and  z  from  m  and  n  alone.  We  therefore  make  the  same  assumption 
generally  made  for  correction  of  large-scale  geometric  distortions  in  sonographs:  that  the 
variations  in  seabed  elevation  can  be  disregarded,  which  corresponds  to  assuming  a  flat, 
horizontal  bottom.  In  that  case  the  variable  z  can  be  eliminated,  and  the  geometric 
mapping  is  completely  described  by  the  two  functions, 

x  =  x  j[m,n]  and  y  =  y*[m,n], 

where  x,[m,n]  and  y,[m,n]  are  the  coordinates  of  the  point  on  the  seabed  that  was 
sampled  to  produce  the  pixel  located  at  [m,n]  in  the  image.  Another  simplification  is 
to  assume  that  the  image  intensity  corresponds  exactly  to  the  backscattering  coefficient, 
since  we  are  not  concerned  here  with  the  intensity  transformation  itself.  We  can  then 
write 

s[m,n]  =  b(x,y)  .  (2.1) 

x  =  x*[m,n] 

y  =  y.[m,n] 

For  the  sonograph  to  provide  an  undistorted  representation  of  the  bottom,  xa[m,n] 
and  ya[m,n]  would  have  to  be  linear  functions  of  m  and  n  with  identical  linear  coeffi¬ 
cients.  In  practice,  however,  they  turn  out  to  be  nonlinear,  due  to  both  the  geometry  of 
the  physical  system  and  the  effect  of  towfish  instabilities.  In  the  next  section,  explicit 
formulas  will  be  derived  for  xa[m,  n]  and  y«[m,n],  involving  parameters  describing  the 
position  and  orientation  of  the  towfish. 

Since  xa[m,n]  and  y*[m,n]  are,  in  general,  nonlinear  functions,  Eq.  (2.1)  can  be 
interpreted  as  describing  a  non-uniform  sampling  of  6(x,y).  In  fact,  plotting  the  sample 
locations  on  the  seabed  plane  does  not  result  in  a  perfect  rectangular  grid,  but  rather  in  a 
sampling  pattern  such  as  the  one  shown  in  Fig.  2.2(a).  This  figure  depicts  the  case  where 
the  towfish  follows  an  ideal  straight  trajectory  at  a  constant  speed  and  with  no  pitching 
or  yawing.  The  dashed  line  running  vertically  in  the  figure  represents  the  projection  of 
the  towfish  trajectory  on  the  bottom,  usually  called  the  bottom  track.  The  other  dashed 
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I  a 

Figure  2.2:  Sampling  patterns  on  the  bottom,  in  the  ideal  case  (a)  and  in  the  presence 
of  towfish  instabilities  (b). 

lines  running  perpendicularly  to  the  bottom  track  represent  the  intersection  of  the  axial 
plane  of  the  transmitting  beams  with  the  seafloor.  These  have  fixed  positions  in  space 
because  each  sound  pulse  is  emitted  at  a  certain  instant  in  time  and  then  propagates 
along  the  direction  that  the  sonar  beam  had  at  the  instant  the  pulse  was  fired.  On  the 
other  hand,  the  return  from  each  pulse  is  received  over  a  period  of  time  during  which 
the  receiving  beam  is  not  stationary  but  rather  keeps  moving  along  with  the  towfish  at 
a  constant  speed.  The  effective  beam  pattern  is  the  product  of  the  transmitting  and 
the  receiving  beam  patterns,  and  if  these  are  the  same,  which  is  usually  the  case  for 
side-scan  sonars,  then  the  effective  beam  can  be  assumed  to  lie  on  the  bisecting  plane  of 
the  axial  planes  of  the  transmitting  and  receiving  beams.  In  other  words,  each  point  on 
the  axial  plane  of  the  effective  beam  is  assumed  to  be  equidistant  to  the  axial  planes  of 
the  transmitting  beam  (which  is  fixed  for  each  pulse)  and  of  the  receiving  beam  (which 
keeps  moving  along  with  the  towfish).  As  a  result,  during  the  time  between  successive 
sonar  pulses  the  effective  beam  scans  the  bottom  at  half  the  speed  of  the  towfish.  Each 
time  a  pulse  is  emitted,  however,  the  effective  beam  jumps  ahead  to  the  location  of  the 
new  transmitting  beam,  so  that  its  average  speed  still  equals  that  of  the  towfish. 

As  the  effective  beam  moves  across  the  seabed,  the  signal  returned  from  the  bottom 
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is  sampled  at  a  constant  rate  as  it  reaches  the  transducers.  The  resulting  samples  come 
from  points  on  the  bottom  that  are  equally  spaced  in  the  along-track  (y)  direction,  since 
the  effective  beam  is  assumed  to  move  at  a  constant  speed  between  pulses.  However, 
the  spacing  between  sampling  points  in  the  cross-track  (x)  direction  varies  as  a  result 
of  the  slant-range  effect.  As  a  consequence,  the  sampling  points  lie  on  a  path  that  is 
initially  curved,  but  approaches  a  straight  line  at  longer  ranges,  with  the  decrease  of  the 
slant-range  effect,  as  shown  in  Fig.  2.2(a).  Part  (b)  of  the  figure  shows  what  the  sampling 
pattern  on  the  bottom  might  look  like  when  the  towfish  is  subject  to  translational  and 
rotational  instabilities.  Because  the  bottom  is  not  sampled  on  an  ideal  rectangular  grid 
but  on  an  irregular  pattern,  the  resulting  image  will  suffer  from  geometric  distortions, 
in  addition  to  those  caused  by  the  slant-range  effect.  The  sampling  interval  along  each 
scan  line  and  between  corresponding  points  of  successive  scan  lines  varies  from  point  to 
point.  Areas  where  the  sampling  interval  is  smaller  will  look  “stretched”  with  respect 
to  areas  where  the  sampling  interval  is  larger.  Furthermore,  the  image  will  look  skewed 
whenever  scan  lines  are  shifted  in  the  cross-track  direction.  The  figure  also  indicates  how 
the  curves  where  the  sampling  points  are  located  may  cross,  resulting  in  backscanning. 


2.2  Mathematical  Model  for  the  Geometric  Distor¬ 
tions 

In  this  section  two  models  will  be  derived  from  the  geometry  of  the  problem.  The 
first  one  will  express  the  absolute  position  of  the  sampling  points  on  the  bottom, 
(x,[m,n],y,[m,n]),  as  a  function  of  the  instantaneous  values  of  the  towfish  attitude 
parameters,  and  will  be  derived  with  respect  to  the  fixed  coordinate  system  (x,y,x).  If 
measured  values  of  the  attitude  parameters  were  available  from  sensors  mounted  on  the 
towfish,  we  would  need  only  substitute  these  measurements  into  this  first  set  of  equa¬ 
tions  in  order  to  calculate  the  coordinates  of  the  sampling  points  on  the  bottom  and  then 
correct  the  geometric  distortions  by  resampling  the  image.  In  our  case,  such  measure¬ 
ments  of  the  attitude  parameters  are  not  available  and  we  will  have  to  seek  an  alternative 
way  of  estimating  the  sampling  point  coordinates.  To  that  end,  a  second  model  will  be 
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Figure  2.3:  The  position  of  the  towfish  is  denoted  by  (x/,y/,z/),  its  pitch  angle  by  <t> , 
and  its  yaw  angle  by  0. 

derived  expressing  the  relative  position  of  the  sampling  points  as  a  function  of  a  set 
of  parameters  which  are  closely  related  to  the  increments  of  the  attitude  parameters 
between  successive  firings.  This  second  model  will  be  used  in  Chapter  4  for  estimating 
the  attitude  parameters. 

2.2.1  Location  of  Sampling  Points 

The  geometry  of  the  problem  is  depicted  in  Fig.  2.3.  The  origin  of  the  coordinate 
system  (x,  y,  z)  is  located  on  the  seabed  plane,  with  y  pointing  in  the  average  direction 
of  the  trajectory  of  the  towfish,  x  pointing  to  the  starboard  side,  and  z  in  the  vertical 
direction.2  The  position  of  the  towfish  at  time  f  is  (x/(f),y/(f),x/(f))  and  its  orientation 
is  described  by  its  pitch  and  yaw  angles,  denoted  ^(f)  and  0(f),  respectively.  The  pitch 
angle  4>(t)  corresponds  to  the  elevation  measured  with  respect  to  the  horizontal  seabed 
plane  and  the  yaw  angle  0(f)  corresponds  to  the  azimuth  measured  counterclockwise 
from  the  y  direction.  It  will  not  be  necessary  to  consider  the  towfish  roll  angle,  since,  as 
explained  in  Chapter  1,  rolling  does  not  contribute  to  geometric  distortions.  Time  f  =  0 
corresponds  to  the  instant  at  which  acquisition  of  the  sonograph  was  initiated,  and  thus 

’Although  in  principle  this  coordinate  system  may  have  any  arbitrary  orientation,  it  is  useful  to  align 
the  y  direction  with  the  trajectory  of  the  towfish  so  that  it  may  be  referred  to  as  the  along-track  direction 
and  x  as  the  cross-track  direction. 
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Figure  2.4:  Relation  between  a  sampling  point  on  the  seabed  and  the  location  and 
orientation  of  the  effective  beam. 


(x/(0),y/(0),z/(0))  is  the  starting  point  of  the  towfish  trajectory. 

Assume  a  sonar  pulse  is  transmitted  at  time  t  =  and  that  the  returned  signal  is 
sampled  at  time  t  =  tR.  The  effective  beam  at  time  tR  is  the  product  of  the  transmitting 
beam  at  time  tj  and  the  receiving  beam  at  time  </*,  and,  as  a  consequence,  its  axial  plane 
is  located  between  the  axial  planes  of  the  transmitting  and  receiving  beams.  Therefore, 
the  orientation  of  the  axial  plane  of  the  effective  beam  at  time  tR  is  completely  defined 
by  the  angles 


Mtn)  4  + 


„d  e.(tR)  4  £(M±*M 


(2.2) 


which  we  call  the  effective  pitch  and  yaw  angles,  respectively.  Additionally,  that  plane 
contains  the  point  (xe(<*),ye(<R),ze(<fl))»  where 

*.(<«)  4  z-iM±±iM'  4  +J(/M  Md  ,.((„)  4 


(2.3) 

In  Fig.  2.4  the  effective  beam  is  represented  by  a  triangular  section  of  its  axial  plane. 
For  simplicity,  only  the  starboard  beam  is  shown.  The  point  of  the  seabed  from  which 
the  signal  sampled  at  time  tR  was  reflected  is  the  sampling  point  (x,(<«),y,(f/i)).  For 
now  x ,  and  yt  will  be  considered  as  continuous  functions  of  time.  This  sampling  point 
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is  located  on  the  line  defined  by  the  intersection  of  the  axial  plane  of  the  effective  beam 
with  the  seabed  plane.  The  projection  of  point  (xe(t/j),ye(t/i),  ze(tfi))  onto  that  line  is 
denoted  by  ( x0(tp),y0(tfi )),  and  the  distance  between  this  point  and  the  sampling  point 
(xs(tR),ys(tR))  is  denoted  by  d.  From  the  geometric  arrangement  depicted  in  Fig.  2.4  we 
see  that  the  equations  we  seek  are 


x,  (tR)  =  x0{tR)  +  dcosde(tR)  (2.4a) 

Vs(iR)  =  y0(tR)  + d  sin  $e(tR).  (2.4b) 

Now  we  need  only  write  x0(tR),y0(tR)  and  d  as  a  function  of  the  towfish  attitude  param¬ 
eters  at  times  tj  and  Ir. 

From  Fig.  2.4  it  is  easy  to  see  that 


*o(<fl)  =  xe{tR)  -  zt(tR)  tan  <t>e(tR)  sin  0e(tR) 
Vo(tR )  =  ye(<*)  +  Zt{tR)  tan  <f>e{tR )  cos  0e(tR). 


As  for  d,  in  Appei  d’  B  it  is  shown  that3 

ef  —  gr2  ±  2r^/l6r4  -f-  4(g2  —  e2  —  4 A2)r2  4-  f2  +  2efg  +  4 e2h2 

2(r2  -  e2) 

where 

a  c(tR  -  ir) 

r  —  - 

2 


(2.5a) 

(2.5b) 


(2.6) 


e  =  [xj{tR)  -  xj{tT)\  cos  ee(tR)  -  [y/(tfl)  +  y/(<r)]  sin  06(<r) 


/  =  M<t)  -  x0(<h)]2  +  (y/(<r)  -  y0(</i)]2  +  z){tT) 

-  [£/(**)  -  x0{tR)f  -  (y/(</?)  -  y0(<*)f  -  zj(tR) 

9  =  [*/(<*)  +  xj(tT)  -  2 x0(tR)] cos  9e(tR)  +  [y/(</i)  +  y/(*r)  -  2y0(tR)] sin 0e(tR) 

=  5  (M<t)  -  i.(<b)]2  +  M<r)  -  !/„((«)]’  + 

+  ( Xf(tR )  -  x0(<«)]2  +  [y/(</?)  -  y0(*«)]2  +  2/(<«))  • 

3In  these,  as  in  all  subsequent  equations,  in  any  compound  signs  (i.  e.,  ±  or  ^)  the  upper  sign  applies 
to  the  starboard  side  and  the  lower  one  to  the  port  side  of  the  sonograph. 
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The  sampling  coordinates  corresponding  to  the  points  of  the  digitized  sonograph  can 
be  obtained  from  Eqs.  (2.4a)  and  (2.4b)  by  taking  tr  =  nTj  and  tn  =  nTj  +  mTa ,  for 
n  =  0, 1, . . .  Nn  —  l,m  =  0,  ±1, . . . ,  ±{Nm  —  1),  where  Tj  is  the  sonar  firing  period,  T j  is 
the  sampling  period  used  in  digitizing  the  returned  signals,  Nn  is  the  number  of  lines  in 
the  sonograph,  and  Nm  is  the  number  of  pixels  per  line  on  each  half  of  the  sonograph. 
As  previously  stated,  if  the  towfish  attitude  parameters  are  available,  as  is  the  case  when 
the  towfish  is  equipped  with  attitude  sensors,  the  equations  derived  in  this  section  can 
be  used  to  calculate  the  true  coordinates  of  each  sampling  point  on  the  seabed  and  the 
image  may  then  be  resampled  to  rectify  the  geometric  distortions. 

2.2.2  Sampling  Displacements 

In  the  absence  of  measured  values  of  the  towfish  attitude  parameters,  we  will  have 
to  rely  on  estimates  of  the  local  degree  of  geometric  distortion  obtained  from  the  image 
itself.  The  slant-range  distortion  can  be  easily  corrected,  as  we  will  see  in  Section  5.1, 
and,  therefore,  we  need  only  be  concerned  with  correcting  the  distortions  due  to  towfish 
instabilities.  For  the  remainder  of  this  chapter,  we  will  assume  that  the  slant-range  effect 
has  been  corrected  to  produce  from  the  original  sonograph  a  new  image  s[/,n],  whose 
cross-track  coordinate  is  no  longer  range  to  the  towfish  but  rather  horizontal  distances 
on  the  bottom.  The  I  coordinate  thus  corresponds  to  the  value  in  pixels  of  the  distance 
d  defined  in  Fig.  2.4  and  given  by  Eq.  (2.6),  in  the  same  way  as  the  m  coordinate  of  the 
original  image  corresponded  to  range  to  the  towfish  in  pixels. 

The  coordinates  of  the  sampling  points  on  the  seabed  associated  with  the  pixels  of 
s[l,n]  are  denoted  by  x#[/,n]  and  yf[/,n].  Our  measures  of  geometric  distortion  in  the 
new  image  will  be  estimates  of  the  differences  in  the  values  of  x,[I,  n]  and  yt[l,  n]  between 
sampling  points  at  fixed  values  of  /  from  one  line  of  the  image  to  the  next.  These 
differences  are  denoted  by 

A„x,[I,n]  =  x„[l,n  +  1]  -  x,[/,n]  (2.7a) 

A„y,[/,n]  =  y,[/,n  +  1]  -y,[/,n],  (2.7b) 

and  are  indicated  in  Fig.  2.5.  Henceforth  A„x,[/,  n]  and  Any,[/,  n]  will  be  referred  to 


37 


Figure  2.5:  The  variation  in  sampling  coordinates  from  one  line  of  the  image  to  the  next 
is  described  by  the  sampling  displacements  An*i[J,n]  and  A„y'[/,n],  shown  here  in  the 
(x,y)  coordinate  system. 

as  the  sampling  displacements  in  the  x  and  y  directions,  respectively.  We  will  see  in 
Chapter  3  how  they  can  be  estimated  from  the  image. 

If  the  estimates  of  A [/,  n]  and  Any,[/,  n]  obtained  through  the  procedure  developed 
in  Chapter  3  were  accurate  enough,  then  the  geometrical  analysis  of  the  preceding  section 
would  not  be  necessary  for  the  correction  of  the  geometric  distortions  in  the  sonograph. 
Indeed,  starting  from  an  initial  guess  of  the  value  of  the  sampling  coordinates  for  the 
first  line  of  the  sonograph,  x,[m,  0]  and  y,  [m,  0],  m  =  0,  ±1, . . . ,  ±(Nm  —  1),  the  sampling 
coordinates  of  the  other  lines  could  be  calculated  on  a  line-by-line  basis  by  recursively 
adding  the  estimates  of  sampling  displacements  A„x,[/,n]  and  Any# [/,  n].  In  that  way, 
the  sampling  coordinates  x,[J,  n]  and  y,[/,n]  would  be  calculated  for  all  points  of  the 
sonograph,  and  could  then  be  used  for  correcting  the  geometric  distortions.  However, 
we  will  see  in  Chapter  3  that  the  estimates  of  sampling  displacements  are  not  accurate 
enough  to  be  used  in  this  way,  and  that  in  order  to  achieve  accurate  correction  of  the 
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geometric  distortions  it  will  be  necessary  to  apply  to  these  estimates  a  lower-dimensional 
model  derived  from  the  geometry  of  the  problem.  If,  from  the  estimates  of  sampling 
displacements,  we  then  estimate  the  parameters  of  this  lower-dimensional  model,  the 
decrease  in  the  dimensionality  of  the  model  can  be  expected  to  reduce  the  estimation 
error  of  the  sampling  point  coordinates.  That  amounts  to  applying  to  the  data  our  a 
priori  knowledge  of  the  structure  imposed  on  the  sampling  operation  by  the  side-scan 
sonar  geometry. 

In  order  to  apply  the  model  derived  in  the  previous  section  to  the  measures  of  sampling 
displacements,  it  would  be  necessary  to  combine  Eqs.  (2.7a)  and  (2.7b)  with  Eqs.  (2.4a) 
and  (2.4b),  after  making  the  conversion  from  m  to  l  prompted  by  the  correction  of 
the  slant-range  distortion.  This  operation  would  result  in  a  set  of  non-linear  equations 
expressing  A nxa[l,n]  and  Any^n]  as  a  function  of  the  towfish  attitude  parameters  at 
times  tj  and  £r.  At  that  point,  those  parameters  could  be  assumed  to  vary  linearly 
between  pulse  firings,  so  that  the  sampling  displacements  for  each  line  of  the  sonograph 
may  be  expressed  as  a  function  of  the  increments  of  the  five  attitude  parameters  between 
successive  pulses.  With  that,  the  dimensionality  of  our  model  would  be  reduced  from 
JVm,  the  number  of  samples  per  line  of  the  sonograph,  to  only  five.  Deriving  these  non¬ 
linear  equations,  however,  would  be  quite  involved  and  using  them  for  estimating  the 
increments  of  attitude  parameters  would  pose  a  costly  non-linear  estimation  problem.  A 
more  practical  approach  is  to  derive  an  approximate  linear  model,  whose  parameters  can 
be  more  easily  estimated  through  linear  recursive  estimation  techniques.  Such  a  model 
is  derived  next. 


2.2.3  Linear  Model  for  the  Sampling  Displacements 

The  first  approximation  required  for  deriving  a  linear  model  for  the  sampling  dis¬ 
placements  is  to  assume  that  the  towfish  remains  stationary  during  the  time  it  receives  the 
returned  signals  for  each  pulse.  In  that  case  the  effective  beam  also  remains  fixed  during 
that  interval  of  time  and  all  the  sampling  points  fall  on  the  line  defined  by  the  intersection 
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of  its  axial  plane  with  the  seabed  plane.  We  then  have,  for  nTj  <  tR  <  (n  4-  1)7/, 

<t>{tT )  =  <£(*/?)  =  <t>e{tR )  =  4>[n] 

0(tT)  =  e(tR)  =  ee(tR)  ±  0[n] 

xj(tr)  —  xj(tR)  =  xe(tR)  =  i/[n] 
yj{tr)  =  y  f{tR)  =  ye(tR)  =  y/[n] 
zj(tT)  =  2/(t/j)  =  ze(*n)  =  z/[n]. 

The  point  ( x0(tR),y0(tR ))  also  remains  stationary  for  each  pulse  and  is  now  denoted  by 
(xc[n],  y0[rc]).  With  this,  Eqs.  (2.5a)  and  (2.5b)  become 

x0[n]  =  xj[n\  —  zj[n]  tan<^[n]sin0[rc]  (2.8a) 

J/o[n]  =  yj[n]  +  z/[n]tan<£[n]cos0[n].  (2.8b) 

The  expression  for  the  distance  d  in  this  case  is  considerably  simplified.  In  fact,  by 
making  the  values  of  the  attitude  parameters  at  times  tR  and  tr  coincide,  and  taking 
tR  =  +  mT„  we  see  from  the  equations  defining  the  variables  that  appear  in  (2.6)  that 

we  now  have 

r  =  m(cTJ  2),  e  =  /  =  g  =  0 

and 

h 2  =  (*/[»*]  -  *o[n])2  +  (y/M  -  y0[n])2  +  z)[n}. 

Thus,  h  is  the  distance  from  point  (*/[n],  y/[n],  2/[n])  to  point  (x0[n],y0[nj),  and  corre¬ 
sponds  to  the  distance  the  sonar  pulse  travels  before  it  first  reaches  the  bottom.  This 
distance  can  be  measured  directly  in  the  sonograph  as  the  height  of  the  water  column  in 
each  line.  Henceforth,  we  will  denote  it  by  h[n]  in  order  to  explicitly  indicate  the  time 
dependence.  Proceeding  with  our  analysis,  we  see  that  Eq.  (2.6)  now  reduces  to 

d  =  ±yJm*(cT,/2y  -  h*[n],  (2.9) 

which  is  exactly  the  relation  prompted  by  the  slant-range  distortion.  Finally,  the  expres¬ 
sions  for  the  sampling  coordinates  become 

x,[m,n]  =  x0[n]  ±  yfm2(cT,/ 2)2  -  h2[n]cos0[n]  (2.10a) 

y4[m,n]  =  y0[n]  ±  y/m2(cTtf 2)2  -  h2[n]sin0[n].  (2.10b) 
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The  equivalent  expressions  for  the  sampling  points  of  the  slant-range  corrected  image 


are 


xa[l,  n]  =  x0[n]  +  l(cTa/2)cos  6[n]  (2.11a) 

y,[l,n]  =  y0{n]  + /(cr,/2)sin%].  (2.11b) 

The  first-order  Taylor-series  expansion  of  these  equations  yields  expressions  for  the 
sampling  displacements  as 

Anx,[/,  n]  ss  (x0[n  +  1]  -  x0[n])  —  l(cTa/ 2)  sin  0[n]  (0[n  -f  1]  -  0[n]) 

A„ya[/,n]  %  (y0[n  +  1]  -  y0[n])  +  f(cr,/2)  cos  0[n](0[n  +  1]  -  0[n]) 

or 

A„x,[/,n]  «  Ax0[n]  —  l(cTt/2)  sin  0[n]A0[n]  (2.12a) 

Any,[/,n]  w  Ay0[n]  -f  l(cTa/2)  cos  0[n]A0[n],  (2.12b) 

where 


Ax0[n]  =  x0\n  +  1]  -  x<>[n] 

^y0[n]  =  y„[n  +  1]  -  y0[n] 

A^[n]  =  %  +  l]-%]. 

Equations  (2.12a)  and  (2.12b)  may  be  modified  to  eliminate  the  dependence  on  the 
instantaneous  value  of  the  yaw  angle.  Consider  the  geometrical  arrangement  for  the  nth 
sonar  pulse,  shown  in  Fig.  2.6.  Let  us  introduce  an  auxiliary  rectangular  coordinate 
system,  ( x',y '),  shown  in  the  figure,  with  origin  at  point  (x0[n],y0[n]),  with  the  x'  axis 
located  on  the  line  defined  by  the  intersection  of  the  axial  plane  of  the  effective  beam  with 
the  seabed  plane,  and  with  the  y'  axis  pointing  in  the  direction  of  the  present  heading 
of  the  towfish.  Thus,  the  axes  of  this  new  coordinate  system  are  rotated  with  respect  to 
the  original  coordinate  system  by  an  angle  equal  to  0[n],  the  current  yaw  angle  of  the 
towfish. 
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Figure  2.6:  Sampling  displacements  in  the  ( x ',  y')  coordinate  system,  defined  to  have  the 
same  orientation  as  the  effective  beam  for  each  pulse. 

The  new  and  the  original  coordinate  systems  are  related  through 

x'  =  (x  -  x0(n])  cos  0[n]  +  (y  -  y<>[n])  sin  (2.13a) 

y'  =  iv  ~  VeW)  cos  0[n]  -  (x  -  x„[n])  sin  0[n].  (2.13b) 

Therefore,  the  sampling  displacements  may  be  converted  from  one  coordinate  system  to 
the  other  through 

Anxi{/,  n]  =  A„x#[/,  nj  cos  0[n]  +  A„y,[f,  n]  sin0[n] 

A„yJ[/,  n]  =  Any,[f,  n]  cos  0[n]  -  Anx,[/,  n]  sin  0[n]. 

From  these  equations  and  (2.12a)  and  (2.12b)  we  obtain  the  expressions  for  the  sampling 
displacements  in  the  (x',  y')  coordinate  system, 

A„x'[/,n]  =  Ax'[n]  (2.14a) 

A„yJ[/,n]  =  Ay'[n]  + /(c7;/2)A0[n],  (2.14b) 
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where  Ai'[n]  and  Ay'[n]  correspond  to  Ax0[n]  and  Ay0[rc]  converted  to  the  new  system. 
This  is  the  lower-order  model  mentioned  before.  We  refer  to  Ax'  [n],  Ay'[n]  and  A0[n] 
as  the  distortion  parameters  to  differentiate  them  from  the  towfish  attitude  parameters 
xj,yj,zj,0  and  4>* 

With  the  linear  model  of  Eqs.  (2.14a)  and  (2.14b),  our  task  is  now  reduced  to  the 
estimation  of  three  parameters,  namely  Ax'fn],  Ay'[n],  and  A0[n],  for  each  line  of  the 
sonograph  from  the  estimates  of  sampling  displacements.  We  will  see  in  Chapter  4  how 
that  can  be  accomplished.  The  resulting  estimates  of  distortion  parameters  will  then  be 
used  for  estimating  the  sampling  point  coordinates  for  all  points  in  the  image,  as  explained 
in  Section  5.2.  Figure  2.7  illustrates  the  overall  process  of  estimating  the  sampling 
point  coordinates.  First,  the  height  of  the  water  column  is  measured  in  the  original 
sonograph  s[m,n]  and  is  then  used  for  correcting  the  slant-range  distortion  to  produce 
a  new  image  s[/,  n].  From  this  new  image,  estimates  of  the  sampling  displacements 
A„x'[/,  n]  and  A„y'[/,n]  are  obtained  by  techniques  introduced  in  Chapter  3.  In  order  to 
reduce  the  estimation  error,  the  linear  model  of  Eqs.  (2.14a)  and  (2.14b)  is  then  applied  to 
the  estimates  of  sampling  displacements,  yielding  estimates  of  the  distortion  parameters 
AxJ,[n],  Ay'[n]  and  A0[n].  The  process  is  then  reversed  to  yield  new  estimates  of  the 
sampling  displacements,  which  are  then  converted  to  the  (x,y)  coordinate  system  and 
added  recursively  to  yield  the  final  estimates  of  sampling  point  coordinates.  As  indicated 
in  the  figure,  it  is  also  possible  to  use  the  estimates  of  Ax' [n],  Ay'[n]  and  A0[n]  for 
obtaining  estimates  of  the  towfish  attitude  parameters.  Though  this  is  not  necessary  for 
correcting  the  geometric  distortions,  it  is  convenient  to  obtain  estimates  of  parameters 
that  describe  more  directly  the  location  and  orientation  of  the  towfish.  The  estimation 
of  attitude  parameters  is  discussed  in  Section  4.2.2. 

In  the  next  chapter  we  will  consider  how  to  obtain  the  estimates  of  sampling  displace¬ 
ments  from  which  the  distortion  parameters  and  attitude  parameters  will  be  estimated 
in  Chapter  4. 


*It  might  seem  more  appropriate  to  call  them  simply  the  model  parameters,  but  we  will  save  that 
designation  for  the  parameters  of  the  state  space  model  introduced  in  Chapter  4. 
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Original  sonograph 
s[m,n] 


*/[«]»»/[«]» 


Figure  2.7:  Overview  of  the  technique  for  estimating  sampling  point  coordinates. 
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Chapter  3 


Measures  of  Geometric  Distortion 


in  Sonographs 


In  order  to  estimate  variations  in  the  towfish  attitude  parameters,  it  is  necessary  to 
extract  from  a  sonograph  some  measure  of  local  geometric  distortion.  In  this  chapter 
we  consider  how  to  obtain  such  measures  and  how  to  estimate  from  them  the  sampling 
displacements  A„x'[l, n]  and  Any' [/, n]  introduced  in  Chapter  2,  from  which  in  turn  the 
distortion  parameters  Ax' [n],  Ay'[n]  and  A0[n]  will  be  estimated  in  Chapter  4.  Because 
of  the  way  the  seabed  is  sampled  on  a  raster-scan  basis,  the  predominant  type  of  geometric 
distortion  in  the  cross-track  direction  is  “skewing,”  that  is,  lateral  shifting  of  the  lines 
so  that  they  are  not  properly  aligned.  On  the  other  hand,  in  the  along-track  direction 
the  predominant  type  of  geometric  distortion  is  “compression”  or  “stretching”  of  the 
image  due  to  variations  in  the  spacing  between  scan  lines  on  the  bottom.  Both  types 
of  distortions  will  be  measured  by  cross-correlating  segments  of  adjacent  lines  of  the 
sonograph. 
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3.1  Assumptions  on  the  Statistical  Properties  of 
the  Sonograph  Image 

Estimating  geometric  distortions  in  an  image  from  the  image  alone  requires  that 
certain  assumptions  be  made  about  the  imaged  scene.  The  basic  assumption  that  al¬ 
lows  us  to  obtain  measures  of  geometric  distortions  in  sonographs  is  that  6(x,y),  the 
two-dimensional  backseat tering  function  of  the  bottom  introduced  in  Section  2.1,  is  a 
wide-sense  stationary  random  process  with  isotropic  autocorrelation  function.  Thus,  if 
Rb{ Ax,  Ay)  denotes  the  autocorrelation  function  of  b(x,  y),  we  have 

M Ax,  Ay)  =  Rb  (7(  Ax)2  +  ( Ay)*)  . 

The  assumption  that  the  backscattering  function  is  wide-sense  stationary  may  not 
hold  if  there  are  significant  variations  in  the  morphology  of  the  bottom.  However,  the 
techniques  developed  here  may  still  be  applied  to  such  images  if  they  are  divided  into 
areas  where  the  morphology  is  uniform  enough  to  justify  the  assumption  of  wide-sense 
stationarity.  The  assumption  of  isotropism  is  valid  if  there  is  no  systematic  orientation  of 
the  features  on  the  bottom  in  some  direction.  Thus,  it  may  not  hold  in  areas  containing 
sand  dunes  or  other  formations  shaped  by  undercurrents. 

Assuming  that  the  backscattering  function  is  wide-sense  stationary  and  that  its  au¬ 
tocorrelation  function  is  isotropic,  the  degree  of  geometric  distortion  in  the  image  can  be 
estimated  by  measuring  variations  in  the  shape  of  the  sample  autocorrelation  sequence 
of  small  areas  of  the  image.  As  explained  in  the  next  sections,  the  sampling  displacement 
A„x'[/,n]  is  estimated  from  the  degree  of  skewing  in  the  l  direction,  of  the  local  sample 
autocorrelation  sequence  at  point  [/,n],  and  Any'[/,  n]  is  estimated  from  the  width  of 
the  main  peak  of  that  sequence  in  the  n  direction.  The  sampling  displacements  may  be 
expected  to  be  comparable  to  the  inter-pixel  distance,  and,  consequently,  computing  the 
autocorrelation  sequence  only  at  lags  equal  to  an  integer  number  of  pixels  would  result 
in  too  coarse  a  quantization  of  their  estimates.  To  achieve  greater  accuracy  in  the  es¬ 
timates,  the  image  is  upsampled  in  both  directions  so  that  the  autocorrelation  function 
may  be  computed  at  non-integer  lags.  The  upsampling  operation  is  carried  out  through 
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zero-padding  of  the  discrete  Fourier  transform  of  the  image.  Denoting  by  su[/',  n']  the  up- 
sampled  image,  by  Nu  the  upsampling  factor,  and  by  ^NltN3  the  (Ni  x  A^-point  Fourier 
transform,  we  have1 


—  ^V*JV,,JV„Nb  {^N,,Sn  {5[/,  n]}} 


The  local  autocorrelation  sequence  of  the  upsampled  image  at  point  [/',  n']  is  estimated 
through  the  sample  autocorrelation  sequence  of  a  block  of  pixels  centered  at  [/',  n'].  In 
principle,  the  averaging  could  be  carried  out  over  blocks  extending  over  two  or  more 
lines  of  the  image,  but  single-line  segments  are  used  because  geometric  distortions  can 
be  expected  to  cause  greater  variation  of  the  autocorrelation  in  the  along-track  direc¬ 
tion  than  in  the  cross-track  direction,  as  a  result  of  the  side-scan  sonar  geometry.  To 
compensate  for  local  variations  in  the  mean  value  and  standard  deviation  of  the  image, 
the  sample  correlation  coefficient  between  line  segments  is  used  instead  of  their  inner 
product.  Assume  that  the  line  segment  used  in  computing  the  correlation  coefficient 
extends  for  L  pixels  on  each  side  of  the  point  being  considered  (in  the  original  image, 
before  upsampling).  Then  the  local  sample  mean  and  standard  deviation  of  the  image 
over  that  segment  are  given  by2 

1  k=L 

M.M  =  57TT  r  *»(*'  +  **«,„']  (3-1) 

iL  +  1  ks-L 


and 


TrrrT  E  («.[('  +WV„,nT  ni)2 

+  1  k=-L 


1 

2 


(3.2) 


The  normalized  sample  cross-correlation  between  a  line  segment  centered  at  point  [/',  n'] 


Mn  practice,  only  a  few  rows  and  columns  of  the  image  are  upsampled  at  a  time,  to  reduce  the 
computer  memory  usage. 

3Because  of  the  upsampling  operation,  the  sample  mean  and  standard  deviation  do  not  change  much 
whether  the  sums  are  carried  out  over  a  continuous  line  segment  of  the  upsampled  image  or  only  over  a 
smaller  number  of  pixels  of  that  segment  separated  by  the  upsampling  interval  Nu,  as  in  these  equations, 
resulting  in  computational  savings. 


47 


Figure  3.1:  Calculation  of  the  sample  cross-correlation  sequence  between  line  segments, 
that  is  used  as  an  estimate  of  the  local  autocorrelation  sequence  of  the  image.  Only  every 
Aruth  pixel  is  taken  in  the  segments  indicated  by  braces,  so  that  the  averaging  is  done 
over  (2 L  +  1)  points,  for  computational  savings  (see  Footnote  2  on  page  47). 


and  another  segment  centered  at  [/'  +  A/',  n'  +  An']  is  given  by 


p[A/',  An';  n'j  -  „']*,(/>  +  A/',  n'  +  An']  +  kN*'  ^ 

(su[/'  +  A/'  +  kNu,  n'  +  An']  -  n,[l'  +  A/',  n'  +  An'])] . 


(3.3) 

This  normalized  sample  cross-correlation  sequence  is  used  as  an  estimate  of  the  auto¬ 
correlation  sequence  of  the  image  at  point  {/',  n']  and  at  lags  A/'  and  An'.  A  graphical 
representation  of  its  calculation  is  given  in  Fig.  3.1. 


3.2  Measuring  Distortion  in  the  Cross-Track  Direc¬ 
tion 

According  to  the  linear  model  model  derived  in  Section  2.2.3,  the  sampling  locations 
corresponding  to  points  of  the  slant-range  corrected  image  s[l ,  n]  are  equally  spaced  along 
straight  lines  on  the  seabed  plane.  As  indicated  in  that  section,  this  is  an  approximation 
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Figure  3.2:  Effect  of  lateral  towfish  displacements  on  the  sampling  pattern.  According 
to  the  linear  model  derived  in  Section  2.2.3,  every  sampling  point  on  each  scan  line  is 
shifted  laterally  by  a  distance  equal  to  the  displacement  of  the  towfish  at  that  line. 

that  assumes  that  the  towfish  remains  stationary  while  receiving  the  returned  signals  of 
each  pulse  and  that  the  increment  in  yaw  angle  between  successive  pulses  is  small.  The 
slant-range  correction  also  assumes  that  the  seabed  is  planar  and  horizontal.  Within 
the  limits  of  these  approximations,  distortions  in  the  /  (cross-track)  direction  are  caused 
mainly  by  lateral  displacement  of  the  towfish,  and  translate  into  a  shift  of  the  sampling 
points  from  one  line  to  the  next  that  was  described  in  Chapter  2  by  the  sampling  dis¬ 
placement  A„x'[/,n].  The  resulting  geometric  distortion  in  that  direction  is  in  the  form 
of  skewing  rather  than  compression  or  stretching  of  the  image,  as  illustrated  in  Fig.  3.2. 
The  sampling  pattern  shown  in  that  figure  assumes  that  the  towfish  bearing  remained 
fixed  in  the  y  direction  as  it  swayed  laterally  from  side  to  side. 

In  the  absence  of  geometric  distortions,  the  autocorrelation  of  the  image  equals  the 
autocorrelation  of  the  backseat tering  function  of  the  seabed,  which  is  assumed  to  be 
isotropic.  In  that  case  its  contour  lines  are  concentric  circles.  However,  when  the  towfish 
is  subject  to  lateral  displacements,  the  autocorrelation  function  becomes  skewed  in  the 
cross-track  direction,  and  its  contour  lines  are  turned  into  ellipses,  as  shown  in  Fig.  3.3(a). 
The  sampling  displacement  Anx'[/,n]  determines  by  how  much  line  (n  -|- 1)  of  the  image 
is  shifted  with  respect  to  line  n,  and,  therefore,  is  reflected  on  the  displacement  of  the 
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Figure  3.3:  Sampling  displacements  in  the  cross-track  direction,  Anx' [/,  n],  are  estimated 
from  the  location  of  the  main  peak  of  the  cross-correlation  sequence  p[Al\  An1;  INU,  nNu], 
at  an  along-track  lag  An'  =  Nu. 

maximum  point  of  the  local  autocorrelation  function  at  a  lag  An  =  1  (or  An'  =  Nu,  in 
the  upsampled  image).  Part  (b)  of  the  figure  represents  a  cross-section  of  p\ A/',  An';  n'] 
at  point  [/',  n']  =  [INU,  nNu\  and  at  an  along-track  lag  An'  =  Nu.  The  maximum  point  of 
this  curve  provides  an  estimate  of  the  sampling  displacement  A„x'[/,n],  The  search  for 
the  maximum  is  limited  to  lags  between  —2  and  2  pixels  (—2 Nu  and  2 Nu  pixels  in  the 
interpolated  image),  since  it  is  unlikely  that  the  lateral  motion  of  the  towfish  between 
successive  pulses  would  be  greater  than  the  corresponding  distance  of  ±40  cm.  Thus 

[l,  n]  =  - arg  max  p\Al\  Nv;  INU,  nNu].  (3.4) 

“-2W.SA/'<3N. 

The  (cTt/2Nu)  factor  is  required  for  the  conversion  from  V  (upsampled  image  coordinates 
in  pixels)  to  x'  (distances  on  the  seafloor  in  meters),  and  the  negative  sign  is  included 
because,  if  line  (n  +  1)  seems  to  be  shifted  in  one  direction  with  respect  to  line  n,  then 
the  sampling  points  must  have  shifted  in  the  opposite  direction  on  the  bottom  from  one 
scan  to  the  next. 

If  the  sonograph  of  Fig.  1.2  is  processed  as  explained  in  Section  5.1  to  eliminate  the 
slant-range  distortion,  we  obtain  the  image  shown  in  Fig.  3.4,  from  which  the  estimates  of 
A„x'  [/,  n]  are  calculated  as  described  above.  The  resulting  estimates  are  shown  in  Fig.  3.5. 
The  cross-correlations  were  calculated  with  I  s  3,  i.  e.,  with  a  window  size  of  seven  pixels, 
and  at  lag  increments  of  1/128  pixel.  The  image  seems  to  contain  primarily  the  noise 


Figure  3.4:  Sonograph  of  Fig.  1.2  after  correction  of  the  slant-range  distortion,  as  ex¬ 
plained  in  Section  5.1. 


Figure  3.5:  Estimates  of  Anx's[l,  n]  for  the  sonograph  of  Fig.  3.4.  The  image  intensity  at 
each  point  is  proportional  to  the  estimate  Anx'[/,n]  at  that  point,  according  to  the  tone 
scale  shown  above. 
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Figure  3.6:  Effect  of  variations  in  speed  and  pitch  angle  (a)  and  yaw  angle  (b)  on  the 
sampling  pattern,  according  to  the  linear  model  derived  in  Section  2.2.3. 

resulting  from  the  estimation  error.  However,  at  a  closer  look,  some  horizontal  patterns 
are  perceived  which  can  be  attributed  to  lateral  displacement  of  the  towfish.  There  are 
also  other  patterns  that  clearly  result  from  linear  objects  in  the  original  sonograph,  such 
as  the  cables  and  scours  on  the  bottom. 


3.3  Measuring  Distortion  in  the  Along- Track  Di¬ 
rection 

In  the  along-track  (n)  direction,  distortions  are  caused  by  variations  in  the  towfish 
speed  and  by  pitching  and  yawing,  which  affect  the  sampling  interval  on  the  bottom  by 
disturbing  the  constant  pace  at  which  the  sonar  beam  would  ideally  move.  Figure  3.6 
illustrates  the  effect  of  variations  of  these  attitude  parameters  on  the  sampling  pattern 
on  the  bottom,  again  under  the  assumptions  discussed  at  the  beginning  of  Section  3.2. 
The  resulting  image  looks  “stretched”  in  areas  where  the  sampling  intervals  in  the  along- 
track  direction  are  smaller  and  “compressed”  in  areas  where  the  sampling  intervals  are 
longer.  A  similar  effect  is  observed  in  the  shape  of  the  local  autocorrelation  sequence  of 
the  image,  as  indicated  in  Fig.  3.7(a).  The  width  of  the  main  peak  of  the  local  sample 
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M  (b) 


Figure  3.7:  Along-track  distortions  cause  local  stretching  or  compressing  of  the  autocor¬ 
relation  sequence  of  the  image  in  the  n  direction  (a).  The  distortions  are  reflected  in  the 
width  of  the  main  peak  of  the  sample  cross-correlation  sequence  p  at  a  cross-track  lag 
A/'  =  0,  and  can  be  measured  through  the  positive-  and  negative-lag  correlation  lengths, 
1+  and  L~  (b). 

cross-correlation  sequence  at  A/'  =  0  reflects  the  distortions  in  the  image  by  becoming 
narrower  or  wider  in  different  areas  of  the  image.  The  width  of  the  mam  peak  can  be 
measured  through  the  correlation  lengths  L+[l,n]  and  L~  [/,  n],  which  are  defined  to  be 
the  positive  and  negative  lags,  respectively,  at  which  the  local  sample  cross-correlation 
sequence  p[ 0,  An';  INU,  n7Vu]  falls  to  a  given  threshold  p0,  divided  by  the  upsampling 
factor,  Nu.  Figure  3.7(b)  illustrates  the  definition  of  the  correlation  lengths.  Thus,  we 
have 

p[0,L+[l,n]N»,lNu,nNu]  =  p0  (3.5a) 

p[0,-L-[l,n]Nu-,lNv,nNu]  =  pa.  (3.5b) 

Figure  3.8  shows  enlargements  of  three  regions  of  the  sonograph  in  Figure  1.2,  each 
one  of  size  31  x  31  pixels,  along  with  the  corresponding  cross-correlation  sequences  at 
the  center  point  and  for  A/'  =  0,  and  the  correlation  lengths  in  each  case.  They  were 

calculated  for  a  window  seven  pixels  wide  ( L  =  3)  and  a  threshold  p0  of  0.5.  The 

first  segment  comes  from  an  area  of  the  image  where  there  is  no  perceptible  geometric 
distortion,  whereas  the  other  two  exhibit  noticeable  stretching  and  compression  in  the 
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along-track  direction. 


3.3.1  Relation  Between  Correlation  Lengths  and  Sampling  Dis¬ 
placements 

In  Section  3.2,  we  saw  that  it  is  straightforward  to  convert  the  location  of  the  main 
peak  of  the  sample  cross-correlation  sequence  into  an  estimate  of  the  sampling  displace¬ 
ment  Anx'[/,  nj.  In  fact,  it  is  only  necessary  to  multiply  the  lag  at  which  the  peak  is 
located  by  the  factor  — ( cTg/2Nu ),  as  seen  in  Fig.  3.3.  In  the  along-track  direction,  a 
complication  arises  because  of  the  fact  that  the  width  of  the  main  peak  of  the  sample 
cross-correlation  function  is  not  directly  related  to  the  sampling  displacement  A„j/'[f,  n] 
that  we  seek  to  estimate,  but  rather  to  the  sampling  interval, 

ft!/,  "1  =  v'(A,l',[(,n])2  +  (A„S:[/,n])^, 

as  illustrated  in  Fig.  3.9.  This  problem  will  be  addressed  later  in  this  section,  but  first  we 
will  consider  how  the  correlation  lengths  L+[l,  n]  and  L~[l,  n]  are  related  to  the  sampling 
interval  Pn[l,  n]. 

It  is  reasonable  to  expect  that  the  correlation  lengths  should  be  approximately  in¬ 
versely  proportional  to  the  sampling  interval  at  each  point  of  the  image.  Specifically,  in 
areas  where  the  sonar  beam  scanned  the  bottom  more  rapidly,  with  a  consequent  increase 
in  sampling  interval,  the  sonograph  image  appears  compressed,  with  a  corresponding  de¬ 
crease  in  correlation  lengths.  Conversely,  in  areas  where  the  sampling  interval  was  smaller 
the  sonograph  appears  stretched  and  the  correlation  lengths  should  be  longer. 

Let  us  consider  for  now  only  the  positive-lags,  for  the  sake  of  simplicity.  Let  Lb[l,n] 
denote  the  correlation  length  of  the  autocorrelation  function  of  the  seabed  backscatter- 
ing  function  b(x,y)  at  point  (x,[/,  n],y,[/,  n]),  defined  in  the  same  manner  as  L+[/,n]  was 
defined  for  the  interpolated  sonograph  image  in  the  preceding  section.  Because  the  au¬ 
tocorrelation  function  of  b(x,y)  is  assumed  to  be  isotropic,  its  correlation  length  equals 
Lb[l,n]  in  all  directions.  L+[l,n]  is  the  distance  in  sonograph  coordinates  that  we  have 
to  travel  from  point  [/,n]  in  the  positive  n  direction  before  the  sample  cross-correlation 
sequence  drops  below  the  threshold  p0  (see  Eq.  (3.5a)).  Similarly,  Lb[l,  n]  is  the  distance 
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Figure  3.8:  Examples  of  the  effect  of  geometric  distortions  on  the  correlation  length,  for 
15  x  15-pixel  areas  where  the  sampling  rate  was  average  (a),  greater  than  average  (b), 
and  less  than  average  (c). 


X 


Figure  3.9:  Relation  between  the  sampling  interval,  Pn[/,  n],  and  the  sampling  displace¬ 
ments,  Anx'[/,  n]  and  Any'[/,  n]. 

the  sonar  beam  has  to  travel  on  the  seabed  from  point  (x,[/,n],y,[/,  n])  before  the  au¬ 
tocorrelation  function  of  6(x,y)  falls  below  that  same  threshold.  Assuming  for  now  that 
I+[/,n]  is  an  integer  number  of  pixels,  one  could  say  that  Lb[l,n\  equals  the  sum  of  the 
next  sampling  intervals  Pn[I,  n]  from  point  (x,[/,n],y,[/,n]),  or 

Lb[l,n]  =  Pn[l,n  +  i ].  (3.6) 

«=o 

In  general,  however,  the  correlation  length  may  not  be  an  integer  because,  as  explained  in 
Section  3.1,  the  correlation  coefficients  are  computed  at  non-integer  lags  so  that  L+[l,n] 
would  not  be  subject  to  coarse  quantization.  Equation  (3.6)  can  be  modified  to  accom¬ 
modate  non-integer  values  of  L+[l,n]  by  taking 

,  .  lUMj-i 

Z,M=  (i+[(.nl- lLi[Z,»U)j>.[J,»+l£;i/,nlJ]+  E  +  (3.7) 

i=0 

where  |Z+[/,  nJJ  denotes  L+[l,n]  rounded  towards  zero. 

Note  that  this  equation  may  also  be  written  as 

£.m 
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where  P+[/,nj  represents  the  average  of  the  £+[/,  n]  sampling  intervals  after  point  [/,  n]. 
This  shows  that  if  Lt[i,n]  is  constant,  or  at  least  fairly  uniform  throughout  the  sono- 
graph,  then,  as  expected,  the  positive-lag  correlation  length  at  point  [/,  n]  is  inversely 
proportional  to  the  sampling  period,  or,  more  accurately,  to  the  average  value  of  the 
sampling  periods  immediately  after  that  point. 

Now,  let  Lno{l ]  represent  the  average  value  of  the  correlation  lengths  in  the  absence 
of  motion  instabilities,  i.  e.,  with  the  to  whs  h  moving  at  a  constant  speed  and  with  the 
heading  aligned  with  the  trajectory.3  No  superscript  is  used  in  Lno[l\  because  it  represents 
the  averages  of  both  X+[/,n]  and  £"[/,  n],  which  should  be  equal  if  there  is  no  geometric 
distortion.  In  the  absence  of  motion  instabilities  the  sampling  periods  equal  vT/,  where 
v  is  the  speed  of  the  deploying  vessel,  and  Eq.  3.7  yields 

Lb[l,n]  =  vTfLno[l]. 

This  allows  us  to  rewrite  that  equation  as 

,  x  lUMJ-1 

vTMH  =  (i+M  -  liJMj)  P.[Z,n  +  IIIMJ]  +  £  />„[/,»  +  i).  (3.8) 

«=0 

This  equation  gives  the  relation  between  the  sampling  interval,  P„[/,n]  and  the 
positive- lag  correlation  length,  L+[l,n\.  Let  us  now  consider  how  the  correlation  lengths 
may  be  related  to  the  sampling  displacements  A„yi[/,n],  the  quantities  we  want  to  esti¬ 
mate.  Recalling  that  Pn[l,  n]  =  \J ( A„xi[/,  n])2  +  ( A„y£  [/,  n])2,  and  since  Anx'  [/,  n]  can  be 
estimated  separately  as  shown  in  the  preceding  section,  we  see  that  Eq.  (3.8)  may  be  used 
for  estimating  A„y'(/,n]  from  the  correlation  lengths  L+[l,n),  provided  the  value  of  Lno[l] 
is  known.  The  estimation  process  would  be  complicated  by  the  fact  that  this  is  a  nonlinear 
equation.  Notice,  however,  that  if  Anx'[/,n]  is  zero,  then  we  have  P„[/,  n]  =  |A„y'[/,n]| 
and  Eq  (3.8)  becomes  a  linear  equation  relating  the  absolute  value  of  A„y'[l,  n]  to  Lno[l). 
This  suggests  that  the  task  of  estimating  AnyJ[/,  n]  can  be  considerably  simplified  if,  be¬ 
fore  calculating  the  correlation  lengths,  we  process  the  image  to  eliminate  the  distortions 

sStrictly  speaking,  the  average  correlation  length  should  be  denoted  Ln(>[nt],  since  the  sonar  beam 
pattern  causes  it  to  vary  with  range  rather  than  with  the  horizontal  distance  on  the  bottom,  as  we  will 
see  in  the  next  section.  However,  we  denote  it  by  Ln0[{]  for  the  sake  of  simplicity,  with  the  conversion 
from  m  to  /  left  implicit. 
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in  the  cross-track  direction,  and  thus  eliminate  the  A„x'[/,n].  That  can  be  accomplished 
by  shifting  the  lines  of  the  image  s[l,n]  by  the  estimated  values  of  the  sampling  dis¬ 
placements  Anx'[/,n],  thus  obtaining  a  new  image  ![/,n]  for  which  A„x'[/,n]  may  be 
assumed  to  be  zero.  This  procedure  will  be  described  in  greater  detail  in  Section  5.1. 
If  the  correlation  lengths  L+[/,n]  are  calculated  as  described  before  from  the  new  image 
s[l,  n],  instead  of  from  s[/,n],  then  A„x'[/,n]  can  be  set  to  zero  in  Eq.  (3.8),  yielding 

vT/L^ll]  =  (ZJM  -  LiitMJ  |Any'[l,n  +  n]J]l  +  E  lAnsilZn  +  i]|, 

1=0 


This  is  the  desired  expression  relating  the  correlation  lengths,  T+[/,  n],  to  the  sampling 
displacements  in  the  along-track  direction,  Any* [/,n].  If,  for  convenience,  we  define 


g[/,n]  = 


lAny^[/,n]| 
vTjLm[1\  ’ 


(3.9) 


then  the  last  equation  may  also  be  written  as 

(Z+[/,n] -[!+[/,  n]j)9[/,n+ LI+[/,n]J]-f  £  q[l,  n  +  i]  -  1.  (3.10a) 

t=0 

Similarly,  for  the  negative-lag  correlation  lengths,  we  obtain 

£,*[/, n]-l 

(^nM  “  L-£nMJ]+  ?[*,"  +  *]  =  L  (3.10b) 

t=0 


At  each  cross-track  distance  /,  (3.10a)  and  (3.10b)  provide  a  set  of  2(Nn  —  1)  indepen¬ 
dent  linear  equations  in  the  ( Nn  —  1)  unknowns,4  g[/,n],n  =  0, 1,2, . . .  ,./Vn  —  2.  These 
equations  can  tk  be  used  to  find  the  deterministic  least-squares  estimates  of  those 
quantities  from  the  correlation  lengths  Z/+[/,n]  and  L~[l,n\  calculated  from  the  image 
![/,  n].  If  for  each  range  index  /  we  define  the  vectors 


«M] 

«['.  i] 


and 


2(JV„  -  1), 


(3.11) 


4There  are  only  (Nn  —  1)  unknowns  because  we  are  dealing  with  sampling  increments  between  lines 
of  the  image,  and  the  number  of  increments  is  one  less  than  the  total  number  of  lines,  N„ . 
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and  form  a  matrix  D[l ]  with  the  coefficients  of  the  q[l,  n]  in  Eqs.  (3.10a)  and  (3.10b), 
then  we  may  write 

mm = m 

The  estimates  of  <?[/,  n],  n  =  0, 1, . . . ,  Nn  —  1  are  finally  obtained  as  the  least-squares 
solution  to  this  equation,  i.  e., 

jM  =  (Hr[WJ)-,/>Til]i!(')-  (3.12) 

This  procedure  amounts  to  inverse  filtering  the  $[/,  n]  along  each  column  of  the  sonograph, 
to  compensate  for  the  time- varying  averaging  implied  by  Eqs.  (3.10a)  and  (3.10b) 

We  now  have  estimates  of  the  quantities  g[/,n]  =  |Any'[/,n]|/uT/Tm)[/].  In  order  to 
obtain  the  desired  estimates  of  the  sampling  displacements  A„y'[/,  n),  it  is  necessary  to 
know  Ln<W],  i-  e.,  what  the  correlation  lengths  should  be  in  the  absence  of  geometric 
distortions.  The  next  section  will  present  a  technique  for  its  estimation. 

3.3.2  Correlation  Lengths  in  the  Absence  of  Geometric  Dis¬ 
tortion 

The  correlation  lengths  are  affected  by  the  composition  of  the  seabed  in  the  area 
scanned  to  form  the  sonograph.  In  fact,  measures  similar  to  the  correlation  lengths 
defined  here  have  been  used  exactly  for  the  purpose  of  classifying  the  different  types  of 
geological  formations  on  the  bottom,  such  as  sand  dunes,  rock  outcrops,  mud,  gravel, 
etc.  [54]  It  is  therefore  necessary  to  know  the  correlation  lengths  associated  with  the 
different  kinds  of  materials  appearing  in  a  sonograph,  if  one  is  to  use  variations  in  the 
correlation  length  as  a  measure  of  geometric  distortion.  This  may  require  some  sort  of 
image  segmentation  process  to  divide  the  image  in  areas  where  the  correlation  lengths 
may  be  considered  uniform.  The  sonograph  we  have  been  considering,  however,  presents 
fairly  uniform  topography  and  does  not  require  such  an  operation. 

The  correlation  lengths  are  also  affected  by  small-scale  variations  caused  by  individual 
objects  on  the  bottom,  such  as  rocks,  scours,  cables,  dunes,  etc.,  and  their  acoustic  shad¬ 
ows.  What  is  important  is  that  underneath  these  small-scale  variations  it  is  possible  to 
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perceive  large-scale  variations  caused  by  geometric  distortions,  which  will  be  used  for  esti¬ 
mating  the  towfish  attitude  parameters.  The  small-scale  variations  in  correlation  lengths 
translate  into  small-scale  errors  in  the  estimates  of  sampling  displacements  A„yJ[/,n], 
which  are  incorporated  as  measurement  noise  in  the  model  derived  in  Chapter  2. 

As  indicated  in  Chapter  1,  the  image  intensity  at  a  certain  point  of  a  sonograph  ideally 
corresponds  to  the  backseat tering  coefficient  of  the  corresponding  point  on  the  seabed. 
Because  of  the  sonar  puke  envelope  and  beam  spread,  however,  the  image  intensity  at  a 
given  point  is  actually  a  weighted  sum  of  the  reflections  from  a  small  area  around  that 
point.  Thus,  a  sonograph  image  in  fact  results  from  the  convolution  of  the  backscattering 
function  of  the  seabed  with  the  sonar  beam  pattern  in  the  along-track  direction  and  with 
the  sonar  pulse  envelope  in  the  cross-track  direction.  If  the  beam  width  changes  as  the 
wavefront  propagates  away  from  the  sonar,  there  are  corresponding  changes  in  the  sample 
cross-correlation  sequence  and,  consequently,  also  in  the  correlation  lengths.  Figure  3.10 
illustrates  the  effect  of  changes  in  beam  width.  The  beam  width  is  approximately  equal 
to  the  length  of  the  transducer  array  in  the  near  field  region  and  then  starts  to  increase 
linearly  with  range  in  the  far  field,  forming  a  wedge-shaped  beam  [11].  On  average,  the 
correlation  lengths  of  the  sonograph  can  be  expected  to  exhibit  a  similar  pattern. 

We  now  need  to  estimate  !„[/]  in  order  to  obtain  the  final  estimates  of  the  sampling 
intervals  Anyi[f,  n].  First  we  select  a  set  AT  of  lines  from  the  sonograph  where  there  is  no 
backscanning  and  where  the  bottom  morphology  appears  to  be  fairly  uniform,  so  that 
//(,[/,  n]  may  be  assumed  constant  and  consequently  the  average  correlation  length  may 
be  a  function  of  /  only.  The  lines  should  also  come  from  areas  of  the  sonograph  where 
visual  inspection  reveals  little  or  no  geometric  distortion,  since  Lno[l]  is  the  correlation 
length  in  the  absence  of  geometric  distortions. 

The  average  of  the  q[l,  n]  over  the  AOv  lines  in  the  set  Af  is* 


n€V 


1  ft][ 

Ntf  ng<V 


5 As  previously  noted,  the  average  correlation  length  is  actually  a  function  of  range  m  rather  than 
horizontal  distance  /.  Therefore,  the  q[l,  n]  should  be  added  not  for  fixed  /  but  for  values  of  /  corresponding 
to  fixed  values  of  m,  according  to  the  slant-range  correction  at  each  line.  To  avoid  making  the  equations 
more  complex,  this  modification  is  again  left  implicit. 
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Figure  3.10:  The  widening  of  the  sonar  beam  pattern  with  increasing  range  is  responsible 
for  a  corresponding  increase  in  the  correlation  lengths  of  the  image. 

JytfVlfLno[l\  n€A/- 

=  v.n  m  £  ( VJ»]  +  l(cr./2)A#W) , 

iSj^vljLnoXl]  „€<v- 

where  in  the  last  step  we  used  Eq.  (2.12).  If  we  symmetrically  add  the  q[/,n]  for 
/  =  ±1,  ±2, . . . ,  ±(JVj  —  1),  the  yaw-dependent  term  disappears,  yielding 

~~  Y  (g[f,  n)  +  q[-l,n])  =  -fi-l-T  m  Y  Ay'[n] 

"W  n€A*  AVvT/Lnol/J  n€Ar 

We  have  thus  reduced  the  problem  to  knowing  the  average  value  of  only  Ay'[n]  (in  the 
( x',y ')  coordinate  system)  over  the  selected  lines.  In  principle,  this  value  is  not  known. 
However,  the  average  value  of  Ay0[n]  (in  the  (ar,  y)  system)  over  the  entire  sonograph 
should  equal  vTj,  the  average  distance  traveled  by  the  towing  vessel  in  the  y  direction 
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between  pulses.  This  assumes  that  at  the  end  of  the  acquisition  of  the  sonograph  the 
distance  between  the  towfish  and  the  deploying  vessel  in  the  along-track  direction  is  the 
same  as  in  the  beginning.  We  will  use  this  assumption  to  calculate  the  average  value  of 
Ay'[n]  over  the  selected  lines  later  on.  For  now  we  represent  this  average  value  by  v'^Tj, 
and  incorporate  it  into  Lno[l]  as  a  scale  factor.  Thus  the  last  equation  yields 


(3.13) 


Figure  3.11(a)  shows  the  average  value  of  q[l ,  n]  for  the  sonograph  we  are  considering, 
calculated  from  51  lines  (out  of  a  total  of  512  lines)  that  appear  to  present  the  least 
amount  of  geometric  distortion.  Since  the  effect  of  the  beam  pattern  may  be  expected 
to  vary  smoothly,  a  fourth-order  polynomial  can  be  fitted  to  that  average  as  seen  in 
the  figure.  That  polynomial  constitutes  our  estimate  of  (v  j  v'j^)Lno\}]-  The  effect  of  the 
beam  pattern  is  seen  clearly  in  this  curve,  added  to  the  constant  level  of  the  correlation 
length  of  the  backscattering  function.  As  expected,  the  curve  is  approximately  constant 
in  the  near  field  portion  of  the  beam  and  then  increases  linearly  with  range  in  the  far 
field,  before  collapsing  as  a  result  of  sound  absorption  by  the  water  as  we  approach  the 
maximum  range  of  100  m.  The  dashed  straight  line  shown  in  the  figure  was  made  to 
coincide  with  the  approximately  linear  portion  of  the  curve.  By  extrapolating  it  to  a 
range  of  zero  we  can  determine  the  value  of  the  correlation  length  of  the  bottom  (times 
the  factor  (v/v#)).  The  resulting  value  is  12.4  cm.  If  we  subtract  that  constant  level 
from  the  polynomial  and  plot  it  twice  “back-to-back”  as  in  part  (b)  of  the  figure,  we 
obtain  a  clear  representation  of  the  effect  of  the  beam  pattern  on  the  correlation  lengths. 

A 

The  estimate  Lno[l ]  may  now  be  multiplied  by  the  estimates  of  q[l,  n]  yielding  estimates 
of  the  scaled  sampling  displacements  (v/uAr)|AnyJ[J,n]|? 


When  we  apply  the  linear  model  of  Eq.  (2.12)  to  these  estimates  of  the  sampling  dis¬ 
placement,  the  resulting  estimates  of  Ay'[n]  and  A0[n]  reflect  the  scaling  factor  (v/v^). 
Then  we  can  finally  calculate  the  value  of  (v/v#)  by  requiring  that  the  average  value  of 
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(b) 


Figure  3.11:  Estimate  of  [v/v^L^H]  obtained  by  averaging  the  correlation  lengths  for 
a  selected  set  of  lines  of  the  sonograph  (a).  The  increment  in  correlation  lengths  due  to 
the  effect  of  the  beam  pattern  is  represented  by  the  curves  in  (b),  which  have  the  shape 
that  could  be  expected  for  the  beam  pattern  of  a  linear  array. 


Ay0[n]  over  the  entire  sonograph  be  equal  to  vTj.  To  keep  the  notation  simple,  we  do 
not  modify  the  equations  of  Chapter  2  to  incorporate  the  scaling  factor  ( v/v '#).  Rather, 
we  ignore  its  existence  for  now  and  simply  rescale  the  parameter  estimates  obtained  in 
Chapter  4  so  that  the  mean  value  of  Ay0[n]  over  the  entire  sonograph  equals  vTj. 

Figure  3.12  illustrates  the  complete  process  of  estimating  the  sampling  displacements, 
and  corresponds  to  one  of  the  blocks  of  the  overall  scheme  shown  in  Fig.  2.7.  Fig.  3.13 
shows  the  estimates  |A„yi[/,n]|  for  our  sample  sonograph,  scaled  by  the  as  yet  unknown 
factor,  (v/v'tf).  The  image  intensity  at  each  point  is  proportional  to  (v/y^)|A„y'[/,  n]| 
according  to  the  tone  scale  shown  below  the  image. 

3.3.3  Detection  of  Backscanning 

The  last  problem  to  address  is  that  we  have  obtained  estimates  not  of  the  sampling 
displacements  A„y'[/,n],  but  of  their  absolute  values.  Therefore,  it  is  necessary  to  detect 
areas  where  A„y'[/,  n]  is  negative,  i.  e.,  areas  where  backscanning  occurs,  and  multiply 
our  estimates  by  minus  one  in  those  areas.  One  possible  way  of  detecting  backscanning 
is  to  look  for  areas  where  objects  appear  in  triplicate.  A  significant  problem  with  that 
approach  is  that  it  is  not  always  easy  to  identify  these  triple  images  because  each  one 
of  them  is  made  from  a  different  viewing  angle  and  may  differ  considerably  from  the 
others.  Fortunately,  there  is  a  more  practical  alternative  based  on  looking  for  indications 
of  backscanning  in  the  estimates  of  A„y'[/,n]  themselves.  If,  for  instance,  the  sampling 
intervals  along  one  line  of  the  sonograph  decrease  with  range  until  they  become  very 
small  and  then  start  to  increase  again,  then  the  portion  of  the  line  after  the  minimum 
was  probably  backscanned  as  a  result  of  yawing.  Visual  inspection  of  the  suspected  area 
can  confirm  whether  backscanning  really  occurred  or  whether  we  have  a  false  alarm. 

Figure  3.14  shows  the  sampling  interval  in  the  along-track  direction  for  a  line  of 
the  sonograph  presenting  backscanning  caused  by  yawing.  Since  the  measured  sampling 
interval  is  typically  very  noisy,  it  was  found  useful  to  fit  a  third-order  polynomial  to  the 
curve  in  order  to  help  determine  the  minimum  point,  sis  shown  in  the  figure.  A  simple 
rule  that  yields  good  results  is  to  decide  that  backscanning  occurred  if  the  minimum 
value  of  the  polynomial  falls  below  a  threshold  chosen  according  to  some  given  criterion. 
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Slant-range  corrected  sonograph 


(v/t>V)A„y;[/,n] 


Figure  3.12:  Overall  technique  for  estimating  sampling  displacements. 
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Figure  3.13:  Estimates  of  (v/v'^Any'gil,  n]|  for  the  sonograph  of  Fig.  3.4.  The  image 
intensity  at  each  point  is  proportional  to  (t>/vjy)|A„y'[/,n]|  at  that  point,  according  to 
the  tone  scale  shown  above. 

For  instance,  the  threshold  may  be  chosen  so  that  the  number  of  false  alarms  and  misses 
is  approximately  the  same.  In  our  case,  a  threshold  of  0.58  pixel  or  approximately  11  cm 
satisfies  that  criterion.  As  for  pitch-induced  backscanning,  if  we  find  two  neighboring 
groups  of  lines  of  the  sonograph  where  the  sampling  interval  is  very  small  on  both  the 
starboard  and  the  port  side,  it  is  possible  that  the  area  between  these  two  groups  of  lines 
was  backscanned  as  a  result  of  pitching.  Again,  visual  inspection  of  the  suspected  strip 
can  indicate  whether  backscanning  actually  occurred. 

Figure  3.15(a)  shows  the  areas  in  the  sonograph  where  backscanning  was  detected,  as 
estimated  through  this  technique.  Visual  inspection  of  the  original  sonograph  reveals  40 
false  alarms  and  40  misses  out  of  the  total  512  lines  (less  than  8  %  in  each  case).  Notice 
that  in  this  sonograph  the  backscanned  areas  alternate  between  the  starboard  and  port 
sides,  indicating  that  the  backscanning  is  due  to  yawing.  Part  (b)  of  the  figure  shows  an 
enlargement  of  the  lower  right-hand  corner  of  the  original  image  for  closer  inspection  of 
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Figure  3.14:  Backscanning  induced  by  yawing  can  be  detected  by  looking  for  lines  where 
the  estimates  |Any'[Z,  n]|  decrease  to  approximately  zero  and  than  start  increasing  again, 
as  shown  in  this  example.  A  third-order  polynomial  is  used  to  help  locate  the  point  where 
backscanning  begins. 

the  backscanned  areas.  The  presence  of  a  cable  on  the  bottom  in  that  area  provides  a 
way  to  verify  the  accuracy  of  the  technique  for  detecting  backscanning.  In  fact,  the  lines 
where  the  cable  appears  to  abruptly  change  direction  were  subject  to  backscanning,  since 
in  all  likelihood  it  ran  smoothly  on  the  bottom  towards  the  lower  right-hand  corner  of  the 
image.  Comparing  the  original  image  with  the  image  where  backscanning  is  indicated, 
we  see  that  the  algorithm  is  very  accurate  in  selecting  the  lines  where  the  cable  appears 
to  change  direction. 

The  last  step  in  the  estimation  of  sampling  displacements  is  to  multiply  |A„y'[/,n]| 
by  -1  in  the  areas  where  backscanning  occurs.  We  thus  obtain  the  estimates  Any'[/,  n], 
which,  along  with  the  estimates  Anx'[/,  n]  obtained  in  Section  3.2,  are  used  in  the  next 
chapter  for  the  estimation  of  the  distortion  parameters  Ax' [nj,  Ay'[n]  and  A0[n]. 
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(b) 


Figure  3.15:  Areas  where  backscanning  was  detected  are  shown  in  black.  In  figure  (a) 
these  areas  are  seen  to  alternate  between  the  starboard  and  port  sides,  indicating  that 
backscanning  in  this  sonograph  was  primarily  caused  by  yawing.  In  (b)  the  lower  right- 
hand  corner  of  the  image  is  enlarged  for  comparison  with  the  original  image. 


68 


Chapter  4 


Estimation  of  Distortion 
Parameters 


In  Chapter  3  we  saw  how  the  sampling  displacements  A„x'[/,n]  and  A„y'[/,n]  can 
be  estimated  from  the  sonograph.  We  now  consider  how  from  those  estimates  we  may 
in  turn  estimate  the  distortion  parameters  Ax'[n],  Ay'[n]  and  A0[n]  of  the  linear  model 
derived  in  Chapter  2.  As  explained  in  Section  2.2.2,  by  applying  a  lower-dimensional 
model  to  the  estimates  of  sampling  displacements,  we  may  expect  a  reduction  in  their 
estimation  error.  As  shown  in  the  overall  scheme  presented  in  Fig.  2.7,  we  will  then  be 
ready  to  estimate  the  sampling  coordinates  of  each  point  of  the  sonograph  and  proceed 
to  rectify  the  geometric  distortions  in  Chapter  5. 


4.1  Selection  of  Observation  Points 

We  now  seek  to  estimate  the  distortion  parameters  Ax'  [n],  Ay'[n]  and  A0[n]  from 
the  estimates  of  A„x'[/,n]  and  Any'[/,r»]  obtained  in  Chapter  3.  The  estimation  is 
carried  out  on  a  set  of  distances  C  =  .  In  principle,  this  set  could  con¬ 

tain  all  columns  of  the  sonograph  from  both  the  starboard  and  port  sides,  i.  e.,  C  = 
{— Nt, . . . ,  —  1, 0, 1, . . . ,  {Ni  —  1)}.  However,  it  has  been  found  advantageous  to  attempt 
to  render  the  estimation  process  more  robust  by  choosing  a  subset  of  the  columns  for 
which  the  estimates  of  sampling  displacements  are  judged  more  reliable  according  to 
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some  criterion.  For  instance,  if  visual  inspection  of  the  sonograph  reveals  areas  where 
the  correlation  lengths  may  vary  markedly  because  of  the  local  morphology,  or  because  of 
artifacts  in  the  image  such  as  acoustic  shadows,  then  the  selected  set  of  distances  should 
be  adjusted,  possibly  on  a  line-by-line  basis,  to  leave  out  those  areas.  We  will  therefore 
represent  the  selected  distances  for  the  nth  line  by  £[n],  As  a  general  rule,  we  have 
found  it  to  be  advisable  to  omit  the  first  few  columns  on  each  side  of  the  bottom  rack 
(usually  about  10  %  of  the  total  number),  since  in  that  area  the  estimates  of  sampling 
displacements  tend  to  be  affected  by  the  undersampling  resulting  from  the  slant-range 
distortion.  That  effect  can  be  clearly  seen  in  the  stripe  running  vertically  in  the  middle 
of  Fig.  3.13,  which  displayed  the  estimates  of  Any'[/,n]. 

Whenever  the  estimates  of  sampling  displacements  reveal  pronounced  yawing  of  the 
towfish,  as  in  the  sonograph  we  have  been  considering,  it  is  also  advisable  to  change  £[n] 
from  line  to  line  so  that  the  selected  columns  always  come  from  the  side  towards  which 
the  towfish  was  turning  at  the  time  that  line  was  acquired.  In  fact,  yawing  causes  the 
beam  to  sweep  the  bottom  faster  on  the  other  side,  which  may  result  in  undersampling 
and  a  consequent  loss  of  accuracy  in  the  estimates  of  sampling  displacements.  Figure  4.1 
illustrates  this  problem  by  presenting  the  estimates  of  A„y'[/,  n]  for  a  line  of  the  sonograph 
in  which  the  towfish  turned  towards  the  starboard  side.  The  straight  line  shown  in  the 
figure  is  the  result  of  fitting  the  linear  model  of  Eq.  (2.14b)  to  the  estimates  of  A„y^[/,n] 
on  the  starboard  side  (/  >  0).  A  good  fit  is  obtained  on  that  side,  but  on  the  port  side 
the  estimates  of  A„y'[/,  n]  depart  from  the  model  because  of  undersampling.  Therefore, 
applying  the  linear  model  to  both  sides  would  result  in  less  accurate  estimates  of  the 
distortion  parameters.  An  estimate  of  which  side  the  towfish  was  turning  towards  for 
each  line  of  the  sonograph  is  easily  obtained  by  computing  the  average  value  of  A„yJ[/,  n] 
on  the  starboard  and  port  sides  of  the  line.  The  side  with  the  smallest  average  sampling 
displacement  is  the  side  the  towfish  was  turning  towards,  and  the  columns  for  £[n]  are 
chosen  from  there. 

For  the  estimations  presented  in  the  next  sections  we  selected  every  fourth  column 
from  columns  number  76  to  475  on  the  side  with  the  smallest  average  sampling  displace- 
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Figure  4.1:  Example  of  a  line  affected  by  undersampling.  The  straight  line  is  obtained 
by  fitting  the  linear  model  to  the  estimates  of  sampling  displacements  on  the  starboard 
side  (7  >  0).  Because  of  undersampling,  the  estimates  of  sampling  displacements  on  the 
port  side  (7  <  0)  are  smaller  than  what  the  linear  model  indicates  they  should  be. 


ment.  Thus 
£[n]  = 


{76, 80,...,  472, 475}, 

{-76, -80,..., -472, -475}, 


475  _  475  

if  £  A„y;[7,n]  <  A„y'[— 7,n], 

/= 76  1=76 

otherwise. 


The  number  of  observations  (p  =  100)  was  limited  by  the  resulting  memory  requirements 
and  execution  times  of  the  algorithms.  In  principle,  better  results  might  be  obtained 
with  a  greater  number  of  observations. 


Figure  4.2  illustrates  the  technique  for  estimating  the  distortion  parameters,  and  is 
part  of  the  overall  scheme  presented  in  Fig.  2.7.  For  each  line  of  the  sonograph  a  set  of  dis¬ 
tances  £[n]  is  first  chosen,  as  described  above.  The  distortion  parameters  for  that  line  are 
then  estimated  from  the  estimates  of  sampling  displacements,  A„x'[7,n]  and  Any'[/,n], 
at  7  €  £[n]  through  one  of  two  techniques  described  in  the  next  sections.  After  all  lines 
of  the  sonograph  have  been  processed,  the  resulting  estimates  of  Ay0[n]  are  then  scaled 
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so  that  their  average  value  equals  one,  thus  compensating  for  the  scaling  factor  ( v/v 
introduced  in  Section  3.3.2.  The  estimates  of  Ax^nj  and  Ay'fn]  are  then  recalculated 

A  A 

to  reflect  the  scaling.  Finally,  the  estimates  of  attitude  parameters,  x/[n],<£[n]  and  0[n], 
may  be  calculated,  though  that  is  not  required  for  the  correction  of  geometric  distortions, 
as  indicated  in  Chapter  2.  These  steps  will  be  examined  in  detail  in  the  next  sections. 


4.2  Deterministic  Least-Squares  Estimation 

One  way  of  estimating  the  distortion  parameters  and  attitude  parameters  is  to  do 
it  separately  for  each  line  of  the  sonograph  using  deterministic  least-squares  estimation. 
This  approach  may  be  used  in  the  absence  of  information  about  the  statistics  of  the 
parameters  and  observation  noise. 


4.2.1  Estimating  the  Distortion  Parameters 

The  model  used  for  estimating  the  distortion  parameters  is 

A^z;[/,n]  =  Ax'[n]  +  vx[/,n]  (4.1a) 

ZQ;[/,n]  =  Ay'[n]  +  /(cr,/2)A0[n]  +  t>y[/,n],  (4.1b) 

which  corresponds  to  Eq.  (2.12a)  and  (2.12b)  with  additional  noise  terms,  vx[/,  n]  and 
uv[/,  n],  that  represent  the  error  incurred  in  replacing  Anx'[/,n]  and  A„yJ[/,  n]  with  their 
estimates. 

Employing  deterministic  least-squares  approximation  amounts  to  choosing  Ax'[n], 
Ay'[n],  and  A0[n]  so  as  to  minimize  the  cost  functions 

=  53  (A^i[/,n]- A^[n])2 

l€£(n] 

e\[n]  =  £  (zQ;[/,n]~  Ay;[n]-/(cr,/2)A^[n])2. 

J€£[n] 

If  p[n]  denotes  the  number  of  observation  points  in  £[n],  and  if  we  define 

Si[n]  =  13  1 

/€£[  n] 
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AnziM,Any;[/,  n] 


Figure  4.2:  Technique  for  estimation  of  the  distortion  and  attitude  parameters. 


3 


Si>M  ^  I? 

l€£[n] 

S*[n]  =  H  A„zi[/,n] 
.'€£[»] 

S„[n]  =  21  £nV',[l,n} 


l€£W 


5lv(n]  =  £  IA„y'[/,n] 

l€£M 

%]  =  p[n]5/a[n]  ~  Sf[n], 
then  the  resulting  estimates  are  given  by 


Ax'[n]  = 

= 

A0[n]  = 


5/a[n]5v[n]  -  5/[n)5/y[n] 
«[n] 

p[nlSit|ti]  -  Si[n]S,[n] 

«!"] 


(4.2a) 

(4.2b) 

(4.2c) 


It  is  now  necessary  to  scale  the  estimates  of  Ax'[n]  and  Ay'[n]  to  compensate  for  the 
{v/v'tf)  factor  introduced  in  Section  3.3.2.  The  first  step  is  to  obtain  estimates  of  the  yaw 
angle  by  cumulatively  adding  the  increment  estimates,  A0[n].  Thus, 

*M  =  0(o]  +  Ea0[4  (4.3) 

i=0 

The  estimate  of  the  initial  yaw  angle,  0[O]  may  be  chosen,  for  instance,  so  that  the 
»  * 
average  value  of  0[n]  over  the  entire  sonograph  be  zero.  However,  the  choice  of  0(0]  does 

not  affect  the  correction  of  geometric  distortions  in  the  image,  but  only  the  orientation 

of  the  corrected  image  as  a  whole.  If  the  assumption  that  the  average  yaw  angle  is  zero 

is  incorrect,  the  only  error  incurred  is  that  the  final  image  will  be  rotated  from  its  true 

orientation  by  an  angle  equal  to  the  difference  between  0[O]  and  the  true  initial  yaw  angle. 

The  next  step  is  to  obtain  estimates  of  x'[n]  and  y'[n]  in  the  (x,  y)  coordinate  system. 

We  first  convert  the  increments  from  the  (x',  y')  to  the  (x,y)  coordinate  system,  using 
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Eqs.  (2.13a)  and  (2.13b),  obtaining 

Ax0[n]  =  Ax'[n]cos0[n]  -  Ay£[n]sin0[n]  (4.4a) 

Ay0[n]  =  Ay'[n]cos0[n]  +  Axj,[n]sin0[n],  (4.4b) 

The  estimates  of  Ay0[n]  must  now  be  scaled  to  make  their  average  value  equal  to  vTj, 
thus  compensating  for  the  (v/x fa)  factor.  In  our  case  t>  =  3  knots  and  Tj  —  0.133  s  (see 
Section  1.5),  resulting  in  vTj  =  20  cm.  Before  scaling,  the  average  value  of  Ay0[n]  is 

22.2  cm,  meaning  that  ( v/v ^-)  equals  1.11.  Equations  (4.4a)  and  (4.4b)  are  now  reversed 
so  that  Ax'[n]  and  Ay'[n]  may  be  recalculated  to  compensate  for  the  {v/v'^)  factor. 
Thus, 

Ax'[n]  =  Ax0[n]  cos  0[n]  +  (t>^/t>)Ayo[n]sin0[n]  (4.5a) 

Ay'[n]  =  (uV/u)^y«[n]cos0[n]  -  Axo[n]sin0[n].  (4.5b) 

Figure  4.3  shows  the  estimates  of  distortion  parameters  obtained  with  this  approach. 
Notice  that  the  estimates  of  Ax'fn]  and  Ay'[n]  have  greater  high-frequency  content  than 
the  estimates  of  A0[n]. 

4.2.2  Estimating  the  Attitude  Parameters 

We  now  proceed  to  obtain  estimates  of  the  towfish  attitude  parameters  from  the 
estimates  of  distortion  parameters.  As  indicated  in  Chapter  2,  this  is  not  required  for 
correcting  the  geometric  distortions  in  the  image,  but  it  is  convenient  to  obtain  estimates 
of  parameters  that  describe  more  directly  the  position  and  orientation  of  the  towfish. 

The  points  (x0[n],y0[n])  are  estimated  by  cumulatively  adding  the  estimates  of  the 
increments  Ax„[n]  and  Ay0[n],  in  the  latter  case  after  multiplication  by  (v'^/v),  i.  e., 

xc[n]  =  xo(0]  +  2^  Ax0[i]  (4.6a) 

tsO 

yoN  =  yo[0]+  £(vV/»)£yo[*J»  (4.6b) 

»=o 

with  the  starting  point  arbitrarily  chosen  as  the  origin  of  the  coordinate  system, 

xo[0]  =  y„[0]  =  0. 
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Figure  4.3:  Estimates  of  distortion  parameters  obtained  through  deterministic  least- 
squares  approximation. 
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To  estimate  the  attitude  parameters,  we  use  Eqs.  (2.8a)  and  (2.8b),  repeated  here  for 
convenience, 


Sofn]  =  x/[n]  —  z/[n]tan<£[n]sin0[n]  (4.7a) 

y0[n]  =  y/[n]  +  z/[n]  tan  ^[n]  cos  0[n].  (4.7b) 

Referring  to  Fig.  2.4,  we  see  that 

zj[n ]  =  h[n)  cos  <f>[n],  (4.8) 

where  /i[n]  is  the  distance  between  the  towfish  and  point  (x0[n],y0[n])  and  can  be  mea¬ 
sured  directly  from  the  sonograph  as  the  height  of  the  water  column  on  each  line.  These 
three  equations  are  all  we  have  to  obtain  estimates  of  X/[n],y/[n],  Zj\n]  and  <f>[n]  from  h[n] 
and  the  estimates  of  x0[n],  y0[n]  and  0[n].  Therefore,  it  is  possible  to  estimate  only  two  of 
the  attitude  parameters,  and  the  third  one  has  to  be  estimated  by  different  means.  One 
solution  to  this  problem  is  to  assume  that  the  along-track  distance  between  the  towfish 
and  the  deploying  vessel  remains  constant  during  the  acquisition  of  the  sonograph.  That 
allows  y/[n]  to  be  estimated  from  navigational  measurements  of  the  deploying  vessel.  If 
such  measurements  are  not  available,  it  is  still  possible  to  resort  to  the  assumption  that 
the  towfish  moved  at  a  constant  speed  v,  equal  to  the  average  speed  of  the  deploying 
vessel.  If  y/[n]  is  measured,  or  assumed  to  satisfy  these  approximations,  Eqs.  (4.7a), 
(4.7b),  and  (4.8)  may  be  solved  for  the  remaining  attitude  parameters,  lateral  position, 
X/[n],  and  pitch  angle,  <£[n].  Thus,  with 


yj\n)  =  nvTj, 

we  obtain  from  Eqs.  (4.7b)  and  (4.8), 

it  i  •  (yo[n)  -  nvT/\ 
4>\n\  =  arcsin  — — - — L  ) . 

\  h[n]  cos  0[n]  / 

Finally,  from  Eq.  (4.7a), 

x/(n]  =  x0[n]  +  fc[n]sin^[n]sin$[n]. 


(4.9) 


(4.10) 


(4.11) 


If  measured  values  of  yj  [n]  are  not  available,  as  in  our  case,  and  the  approximation  of 
(4.9)  is  adopted,  then  variations  in  towfish  speed  from  the  assumed  fixed  value  v  affect 
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the  estimates  of  pitch  angle,  calculated  through  Eq.  (4.10).  Therefore,  this  equation 
should  be  used  only  with  the  clear  understanding  that  the  resulting  estimates  4>{n]  reflect 
not  only  the  towfish  pitch  angle,  but  also  variations  in  its  speed. 

Equations  (4.3),  (4.10),  and  (4.11)  can  be  used  for  estimating 0[n],  <£[n],  and  x/[n]  from 
Ai'[n],  Ay'[n]  and  A 0[n].  These  equations  can  also  be  used  for  determining  approximate 
confidence  intervals  for  the  resulting  estimates.  If  the  estimation  errors,  vx[/,n]  and 
vv[l,  n],  are  assumed  to  be  Gaussian  random  variables,  it  can  be  shown  [5]  that  the 
variances  of  the  estimates  of  distortion  parameters  are  given  by 


Equations  (4.4a)  and  (4.4b)  may  be  used  to  obtain  first-order  estimates  [43]  for  the 
variances  of  Ax0[n]  and  Aye[n], 


=  cos2  %]<r|^(n]  +  sin2  0[n)(T^\n\ 

+  (Ax' [n]  sin  0[n]  +  Ay'[n]  cos  0[n]J  a |[n] 


and 


=  cos2  £[n]<r|-[n]  +  sin2  %]<r|-[n] 

-I-  (Ax' [n]  cos  9[n]  -  Ay;[n]  sin0[n])2  cr)\n). 


Making  the  simplifying  assumption  that  the  increment  estimates  A0[n],  Ax0[n],  and 
Ay0[n]  for  different  lines  of  the  sonograph  are  independent,  we  have 


*|[n]  = 


«=o 
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<(«]  = 

tsO 


n 


By  adding  ±2<Tio[n),  ±2<rie[n],  and  ±2 a6\n\  to  x0[n],y0[n],  and  0[n],  respectively,  we  ob¬ 
tain  confidence  intervals  for  these  estimates.  If  the  assumptions  made  above  are  valid, 
namely,  that  the  estimation  errors  are  Gaussian  random  variables,  that  the  increment 
estimates  for  different  lines  are  independent,  and  that  the  first-order  estimates  of  [n] 
and  <t2-—  [n]  are  accurate,  then  these  confidence  intervals  define  the  ranges  within  which 
x0[n],  yG[n],  and  0[n]  can  be  expected  to  be,  with  95  %  probability. 

By  applying  Eqs.  (4.11)  and  (4.10)  to  the  lower  and  upper  limits  of  the  confidence 
intervals  of  x0[n]  and  y0[n]  we  also  obtain  confidence  intervals  for  X/[n]  and  4>[n\.  The 
resulting  estimates  of  attitude  parameters  and  their  confidence  intervals  are  shown  in 
Fig.  4.4.  In  view  of  all  the  approximations  and  simplifying  assumptions  that  were  made 
in  the  derivation  of  the  confidence  intervals,  they  should  be  considered  only  as  approxi¬ 
mate  indicators  of  the  magnitude  of  the  estimation  errors  of  the  attitude  parameters.  For 
that  reason,  they  are  not  labeled  as  95  %  confidence  intervals  in  the  figures.  However, 
notice  that,  though  the  precise  significance  of  these  confidence  intervals  is  not  rigorously 
established,  they  are  useful  for  indicating  relative  trends  in  the  evolution  of  the  estimates 
and  in  their  relative  accuracy  with  respect  to  one  another.  For  instance,  notice  that  the 
confidence  intervals  grow  larger  with  increasing  n,  since  calculating  the  attitude  param¬ 
eters  involves  cumulatively  adding  the  estimates  of  increments  of  distortion  parameters, 
AxJ,[n],  Ay'[n],  and  A0[n].  Thus,  the  estimation  error  adds  up  and  increases  with  n.  The 
confidence  intervals  also  indicate  that  the  estimates  of  yaw  angle,  0[n],  are  more  accurate 
than  the  estimates  of  pitch  angle,  <^[n].  That  is  explained  by  the  fact  that  during  the 
acquisition  of  our  sonographs,  the  altitude  of  the  towfish  above  the  bottom  was  approx¬ 
imately  one  tenth  of  the  maximum  range  (as  is  commonly  the  case).  Thus,  variations  in 
yaw  angle  cause  a  displacement  of  the  scan  lines  at  long  ranges  that  is  10  times  larger 
than  the  displacement  caused  by  equal  variations  in  pitch  angle.  At  shorter  ranges  the 
difference  is  not  so  marked,  but  still  the  average  effect  of  yawing  over  all  ranges  is  greater 
than  the  average  effect  of  pitching,  resulting  in  better  estimates  of  the  yaw  angle. 
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Figure  4.4:  Estimates  of  attitude  parameters  obtained  through  deterministic  least-squares 
estimation.  The  confidence  intervals  are  indicated  by  dashed  lines. 


4.3  Bayesian  Linear  Least-Squares  Estimation 


Another  approach  to  the  estimation  of  the  distortion  parameters  is  to  use  Bayesian 
least-squares  (or  minimum  variance)  estimation,  which  incorporates  statistical  informa¬ 
tion  about  the  parameters  and  measurement  noise,  to  minimize  the  variance  of  the  es¬ 
timation  error  [43].  Ideally,  the  statistical  properties  of  the  parameter  and  noise  vectors 
should  be  derived  independently,  from  experiments  or  from  a  theoretical  analysis  of  the 
dynamics  of  the  towfish  and  of  the  autocorrelation  of  the  backscattering  function  of  the 
bottom,  respectively.  It  is  also  possible  to  attempt  to  derive  such  information  from  the 
same  data  from  which  the  distortion  parameters  are  estimated,  but  caution  should  be 
exercised  in  the  interpretation  of  the  results  thus  obtained.  This  section  presents  pro¬ 
cedures  for  Bayesian  linear  least-squares  estimation  of  the  distortion  parameters,  first 
from  each  line  separately,  and  then  by  recursively  updating  the  estimates  at  each  line. 
These  procedures  are  described  with  the  intent  of  providing  a  framework  in  which  inde¬ 
pendently  derived  statistical  information  about  the  parameters  and  measurement  noise 
may  be  easily  incorporated  in  the  future. 

4.3.1  Line-by-Line  Estimation 

The  correlation  between  the  distortion  parameters  and  measurement  noise  can  be 
described  more  compactly  in  matrix  form  if  we  write  Eqs.  (4.1a)  and  (4.1b)  as 

£[n]  =  C[n)Qn\  +  j>[n], 
where  £[n],£[n]  and  ji[n]  are  the  vectors 

Ax;[n] 

Ay^nl 

A0[n] 


(4.11) 
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and  C[ri\  is  the  observation  matrix 


1  0  0 

a  10  0 

C[n\  =  .  (4.13) 

0  1  h{cTt/2) 

.  0  1  lp(cT,/ 2)  , 

Let  rri(\n\  and  mjn]  denote  the  expected  values  of  £[n]  and  £[n],  respectively,  and 
assume  that  the  expected  value  of  u[n]  is  zero.  Let  Ajn],A^[n]  and  i?[n]  denote  their 
covariance  matrices.  Further,  let  A^[n]  and  A^v[n]  denote  the  cross-covariance  matrices 
of  £[n]  with  £[n]  and  n[n],  respectively.1  Then  the  linear  minimum  variance  estimate  of 
£[n]  given  £[n]  is  given  by  [43] 

+  A«[n]A^[n]  (C[n]  -  mc[n])  ,  (4.14) 

where,  from  Eq.  (4.11), 

mc[n]  =  C[n]m([n] 

A(t[n)  =  A([n]CT[n|  +  A<„[n] 

A<[n]  =  C!nJA([n]C:r|n!  +  C[r.]A<Jn]  +  Af„[n]CT[nJ  +  B|n]. 

If  numerical  values  for  the  covariance  matrices  are  stipulated  theoretically  or  measured 
from  the  data  itself,  these  equations  may  then  be  used  for  estimating  the  distortion 
parameter  vector  £[n]  on  a  line-by-line  basis. 

4.3.2  Recursive  Estimation 

Bayesian  least-squares  estimation  seeks  to  minimize  the  variance  of  the  estimation 
error  by  exploiting  the  correlation  between  the  components  of  the  parameter  and  noise 
vectors  at  time  n.  However,  it  is  also  possible  to  take  into  account  their  temporal  cor¬ 
relation.  In  other  words,  the  distortion  parameter  vector  £[n]  can  be  estimated  through 

1  In  this  problem,  A^[n]  is  in  general  non-zero  because  the  amplitude  of  the  measurement  noise  can  be 
affected  by  the  degree  of  geometric  distortion,  and  is  therefore  correlated  with  the  distortion  parameters. 
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stochastic  filtering.  In  this  particular  problem,  a  Kalman  filter  is  a  natural  choice,  since 
Kalman  filters  assume  that  the  temporal  correlation  of  the  parameters  to  be  estimated 
can  be  described  by  a  state-space  model,  which  in  our  case  is  a  reasonable  assumption, 
given  that  our  parameters  are  related  to  the  evolution  of  a  hydrodynamic  system.  We 
employ  a  standard  Kalman  filter  to  obtain  a  causal  estimation  procedure,  though  it  is 
also  possible  to  use  smoothing. 

For  an  Nth-order  state  space  model,  the  3-dimensional  vectors  of  distortion  parame¬ 
ters  will  be  stacked  to  form  a  3N-dimensional  state  vector 

£["] 

£[n- 

i[n-(N-  1)].  J 

Adding  to  our  model  an  Nth-order  state  propagation  equation,  we  obtain  the  state-space 
model 

i[n  +  1]  =  Afnjifn]  +  £[n]j£[n]  (4.15) 

£[n]  =  C[n]i[n]  +  n[n],  (4.16) 

where  n>[n]  is  a  white  process  noise  vector.  For  the  second  equation  to  be  equivalent  to 
Eq.  (4.11),  the  new  observation  matrix  is  given  by 

C[n]  =  [C[n]  0  •••  0]. 

Following  the  notation  commonly  adopted  in  the  literature,  the  covariance  matrices  of 
the  noise  vectors  m[n]  and  n[n]  will  be  denoted  by  Q[n]  and  i2[n],  respectively,  and  their 
cross-covariance  matrix  will  be  denoted  by  $[n).  Two  issues  will  be  addressed  next: 
how  to  use  this  state-space  model  to  recursively  estimate  the  parameters,  and  how  to 
determine  the  unknown  matrices  A[n],B[n),  Q[n],  R[n\  and  S[n]  of  the  model. 

One  way  to  determine  the  matrices  of  the  state-space  model  of  Eq.  (4.15)  is  through  a 
theoretical  analysis  of  the  towfish  dynamics,  possibly  in  conjunction  with  the  estimation 
of  some  of  the  parameters  in  the  theoretical  model  from  actual  attitude  measurements 
obtained  through  experiments  in  which  sensors  are  mounted  on  the  towfish.  If  the  model 
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thus  derived  still  contains  unknown  parameters  or  time- varying  parameters  that  depend 
on  the  conditions  during  the  survey  and  cannot  be  determined  beforehand  through  exper¬ 
iments,  then  these  parameters  and  the  state  vector  itself  may  be  concurrently  estimated 
through  adaptive  Kalman  filtering. 

The  general  problem  we  face  is  to  simultaneously  estimate  the  state  vectors  i[n] 
and  a  number  of  unknown  parameters  of  the  state-space  model,  which  will  be  combined 
into  a  parameter  vector  &[«].  In  the  next  sections  we  consider  how  the  state  vectors  and 
model  parameters  can  be  estimated  separately  by  two  recursive  algorithms  and  how  these 
algorithms  can  be  combined  to  form  an  adaptive  Kalman  filter.  Our  treatment  of  these 
algorithms  is  not  exhaustive;  the  intent  here  is  to  provide  a  general  framework  for  the 
application  of  recursive  estimation  and  identification  techniques  to  the  results  derived  in 
Chapters  2  and  3. 


Recursive  Identification  of  the  State  Vector 

The  discrete-time  Kalman  filter  is  the  optimal  solution  to  the  problem  of  recursively 
calculating  the  Bayesian  linear  least-squares  estimate  of  the  present  value  of  a  discrete¬ 
time  stochastic  process  given  the  present  and  past  values  of  another  stochastic  process. 
It  assumes  that  the  dynamics  of  the  process  to  be  estimated  and  its  relationship  with 
the  observation  process  can  be  described  through  a  state-space  model.  The  details  of 
the  derivation  of  the  Kalman  filter  equations  can  be  found  in  [1]  and  [26].  Using  a 
standard  Kalman  filter  would  require  inverting  a  matrix  of  order  2p,  the  dimension  of 
the  observation  vector  £[n],  in  the  measurement  update  step  of  each  iteration.  In  our 
case,  a  is  more  computationally  efficient  to  use  the  so-called  information  filter  version  of 
the  Kalman  filter,  that  requires  the  inversion  of  a  matrix  of  order  3 N,  the  dimension  of 
the  state  vector  £[n].  The  information  filter  uses  the  quantities 

fl[n|n  —  1]  =  f}2nP~x\n\k  —  l]^[n|n  —  1] 

fll»M  =  0*P-l[n\n]*lnH 
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where  P[n|n]  and  P[n|n  —  1]  are  the  variances  of  the  estimation  and  prediction  errors, 
respectively,  i.  e., 


P[n|n]  =  E  |  (i[n|n]  -  =[n])  (“[n|n]  -  i[n])T| 

P{n|n  -  1]  i  p|(i[n|n-l]-i[n])  (i[n|n-l]-i[n])T|, 


and  /?  is  a  forgetting  factor,  0  <  /3  <  1,  used  to  decrease  the  weight  of  past  observations. 
With  /3  =  1  we  obtain  the  standard  Kalman  filter.  In  the  case  of  a  linear  model,  such 
as  the  one  given  in  (4.15),  the  estimates  of  the  state  vector  that  minimize  P[n|n]  can  be 


calculated  recursively  according  to  the  following  equations: 

Measurement  Update  Step 

£[n|n]  =  £[n|n  —  1]  +  Cr[n]/2-1  [n]£[n]  (4.17a) 

P-1[n|n]  =  P-1[n|n  —  1]  +  f3~inCT[n\R~l[n  —  l]C[n]  (4.17b) 

=[n|n]  =  /S~2"P[n|n]£[n|n]  (4.17c) 

Prediction  Step 

P[n]  =  A[n]  —  P[n]5[n]P_I[n]C[n]  (4.18a) 

G[n]  =  )92{n+1)F-T[n]P-1[n|n]P-1[n]  (4.18b) 

H[n]  =  (4.18c) 

fl[n  +  l|n]  =  (/-^[n]PT[n])(p-T[n]a[n|n]  +  G[n]B[n]5[n]p-1[n]^[n])  (4.18d) 

P-1[n  +  l|n]  =  /T2n  (/  -  H[n]BT[nfj  G[n ]  (4.18e) 


These  equations  assume  that  the  process  noise  ia[n]  and  the  measurement  noise  n[n] 
are  white.  If  that  is  not  the  case,  the  state-space  model  can  be  augmented  to  incorporate 
a  model  of  the  temporal  autocorrelation  of  the  noise  processes. 

Recursive  Least-Squares  Identification  of  the  Model  Parameters 

An  autoregressive  (AR)  model  was  adopted  for  the  dynamics  of  the  state  vector. 
While  this  is  not  the  most  general  model  that  could  be  adopted,  it  illustrates  this  tech¬ 
nique  without  unnecessarily  complicating  its  formulation.  In  fact,  an  autoregressive 
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model  can  be  easily  incorporated  into  the  state-space  model  of  (4.15)  by  simply  adopting 
a  companion-form  matrix  ^4[n].  The  TVth-order  autoregressive  model  given  by 


£[n  -I- 1]  =  Aj£[n]  +  A2|jn  -  1]  -I- ...  -I-  As§,[n  -  N  +  1)  +  jyfn], 


corresponds  to  Eq.  (4.15)  with 


A[n\  = 

>li(n]  >42[n]  •••  /!*[«] 

I  0  0  0 

0  /  0  0 

and  i?[n]  = 

1 

0 

...  o 

*~i 

o 

...  o 

_ 1 

0 

where  N  is  the  dimension  of  the  parameter  vector  and  where  I  and  0  represent  identity 
and  zero  matrices  of  dimension  equal  to  the  number  of  parameters  used  in  the  state  space 
model,  in  our  case,  three  (namely,  Ax' [n],  Ay'[n]  and  A0[n]). 

Our  objective  is  to  estimate  the  blocks  A\[n],  .42[n], . . . ,  i4jv[n]  of  matrix  A[n).  It  will 
be  useful  to  combine  the  first  three  rows  of  -4[n],  which  contain  these  sub-matrices,  into 
a  9Ar-dimensional  parameter  vector 


M)rw 
wrw  . 

(A)JW 


(4.19) 


where  (>4),[n]  represents  the  tth  row  of  4[n].  It  will  also  be  useful  to  rearrange  the  state 
vectors  into  a  new  (3  x  9Ar)-dimensional  matrix 

ir[n]  0  0 

0  ir(n]  0  , 

0  0  2r[n] 

which  is  the  multidimensional  counterpart  of  what  is  usually  called  the  regression  vector 
in  the  one-dimensional  case.  The  autoregressive  model  may  be  used  to  predict  the  next 
value  of  the  state  vector  from  the  present  and  (N  —  1)  past  values.  Denoting  by  £[n  + 1  |n] 
the  predicted  value  of  £[n  +  1]  given  f[n],£[n  -  1], . . .  ,£[n  —  N  - 1-1],  we  have 


!("  +  l|n]  =  xMaM- 
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There  are  several  possible  criteria  for  estimating  the  model  parameters  &[n]  from 
the  vectors  £[n]  [38].  One  is  to  minimize  the  mean-square  value  of  the  prediction  error. 
This  criterion  motivates  a  number  of  closely  related  algorithms,  among  which  is  found 
the  weighted  least-squares  method.  In  this  method  the  parameter  estimate  is  chosen 
according  to 

&[n]  =  arg  min  £  02(n-*)  ($[*  +  l]-#  +  l  |*])T  Q~l  [fc]  (i[k  +1]  -  £[k  +  l|/r])  , 

fc=0 

where  /3  is  again  a  forgetting  factor  satisfying  0  <  $  <  1.  The  recursive  version  of  the 
weighted  least-squares  algorithm  consists  of  calculating  a  new  estimate  q[ti]  once  £[n]  is 
available,  and  then  using  it  to  calculate  the  prediction  £[n  +  1  |n]  of  its  next  value.  It  can 
be  shown  that  if  these  estimates  are  to  satisfy  the  least-squares  criterion  given  above, 
they  can  be  computed  recursively  as 

J[n]  =  U[n  -  l]xrN  (02<2[n]  +  x[«]^[«  -  l]xT[n])  1  (4.21a) 

=  fi[n  -  1]  +  J[n]  (^[n]  -  x{n]a[n  -  1])  (4.21b) 

t'N  «  ^  ({/!"  -  1]  -  t'ln  -  llxT  (^QW  +  X[n]^[n  -  (4.21c) 

xWtl[n  -  1])  . 

For  the  derivation  of  the  recursive  least-squares  algorithm  see  [38]. 

The  Extended  Least-Squares  Method 

If  the  model  parameters  are  known,  then  i[n]  can  be  calculated  from  the  observa¬ 
tions  £[n]  through  a  Kalman  filter.  Conversely,  if  the  state  vectors  are  known,  then  a[«] 
can  be  calculated  from  them  and  the  observations  through  the  recursive  least-squares 
method.  Figure  4.5(a)  depicts  this  situation.  When  both  the  state  and  parameter  vec¬ 
tors  are  unknown,  the  Kalman  filter  and  the  recursive  least-squares  algorithm  can  be 
combined  into  what  is  known  as  the  extended  least-squares  method,  as  shown  Fig.  4.5(b). 
It  consists  of  calculating  the  parameter  estimate  &[n]  at  each  iteration  from  the  present 
and  past  state  vector  estimates,  and  then  using  it  in  the  prediction  step  of  the  Kalman 
filter.  The  state  and  parameter  estimates  obtained  through  this  method  are  suboptimal 
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(a) 


£["] 


(b) 

Figure  4.5:  The  Kalman  filter  and  recursive  least-squares  algorithms  (a)  are  combined 
in  the  extended  least-squares  method  (b). 

because  they  are  not  jointly  estimated,  but  they  can  be  shown  to  converge  to  the  optimal 
estimates  under  certain  conditions  [1]. 

Since  the  covariance  matrices  are  unknown,  they  are  estimated  by  recursively  com¬ 
puting  the  average  of  the  outer  products  of  the  residual  errors  from  the  state  and  model 
parameter  estimation  steps.  This  too  is  a  suboptimal  procedure  that,  under  certain  cir¬ 
cumstances,  can  be  shown  to  converge  to  the  optimal  estimates  [40].  The  complete  set 
of  equations  for  the  extended  least-squares  algorithm  is: 
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State  Estimation 

fl[n|n]  =  &[n|n  —  1]  -f  CT[n]i?_1[n  —  l]£[n]  (4.22a) 

P-1[n|nj  =  P-1[n|n  —  1]  +  P~2nCT[n\Brl\n\C\n\  (4.22b) 

»[n|n]  =  /?-2nP[n|n]a[n|n]  (4.22c) 

Model  Estimation 

J[n]  =  U[n  -  l]xT(n|n]  -  1]  +  xN«M«  -  l]xT[«M)  *  (4.23a) 

&[n]  =  a[n  -  1]  +  «7[n]  (|[n|n]  -  x[njn]a[n  ~  1])  (4.23b) 


U[n]  =  jp  (u[n  -  1]  -  U\n  -  l]xT[nM  (/32 Q\n ]  +  x[n|n]B[n  -  l]xT[n|n])  1  (4.23c) 
Xln|n]t/[n  -  1])  . 

Estimation  of  Covariance  Matrices 


&[n] 

£N 

Q\n] 

/R[n] 

5[n] 


|[n|n]-x[n-l|n-l]d[n] 

£["]  -  C[n]i[n|n] 

((4^)  -  ‘I +&M*rM) 

1  _  («-'[" -1]£W)T 

(ti-7)  fS[n ' 11  +  anteT|nl) 


(4.24a) 

(4.24b) 

(4.24c) 

(4.24d) 

(4.24e) 


State  Prediction 


F[n]  =  A[nJ  -  (4.25a) 

G[n]  =  ^n+1>P-r(n]P~1[n|n]F-1(n]  (4.25b) 

B[n]  =  G{n]£[n]  ^BT[n]G[n]B[n] -f  (<2[n]  -  S[n]B_1[n]S[n])  ^  (4.25c) 

fl[n  +  l|n]  =  (/-//(n]BTln])(F-T(n]i[n|n]  +  G[n]P[n]5[n]F'MnK(n])  (4.25d) 

P_1[n  +  ljn]  =  r2n(l-H[n]BT[n])G[n]  (4.25e) 

Figure  4.6  shows  the  estimates  of  distortion  parameters  obtained  through  this  algo- 
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rithm,  with  N  =  4  and  0  =  1.  The  initial  state  estimate  was 

0 

i[0|  -  1]  =  vTj  , 

0 

and  initial  estimates  for  the  covariance  matrices  were  obtained  from  the  parameter  and 
error  estimates  of  the  deterministic  least-squares  algorithm.  The  initial  estimates  of 
the  state-space  model  parameters  corresponding  to  the  autoregressive  coefficients  of 
Ax' [n],  Ay'[n],  and  A0(n]  were  equal  to  1/N  and  all  the  coefficients  for  the  cross  terms 
were  equal  to  zero.  Thus,  recalling  the  definition  of  a[n]  in  Eq.  (4.19),  we  have,  for 
N  =  4, 

r  ,  0.25,  for  1  =  1,4,7,10,14,17,20,23,27,30,33,36, 

<M“1]  = 

0,  otherwise. 

With  this  choice  of  the  dr,[n]  the  initial  autoregressive  model  amounts  to  taking  the 
prediction  of  the  next  value  of  the  state  vector  as  the  average  of  its  last  N  values. 
Finally,  the  initial  covariance  matrix  t/[ — 1]  was  a  9iV-dimensional  identity  matrix. 

The  state-space  model  causes  the  estimates  of  Ax' [n]  and  Ay'[n]  to  have  less  high- 
frequency  content  than  the  corresponding  estimates  obtained  through  the  determinis¬ 
tic  least-squares  algorithm.  The  estimates  of  A0[n],  however,  do  not  change  as  much 
since  they  already  have  low  high-frequency  content  in  the  deterministic  least-squares 
case,  probably  as  a  result  of  the  smaller  estimation  errors  achieved  for  that  parameter. 
Figure  4.7  shows  the  estimates  of  attitude  parameters  for  the  extended  least-squares 
algorithm.  These  do  not  differ  significantly  from  the  estimates  obtained  through  de¬ 
terministic  least-squares  estimation,  which  could  be  expected  since  the  only  a  priori 
information  introduced  was  the  autoregressive  structure  of  the  state-space  model,  and 
its  order.  If  some  of  the  model  parameters  are  determined  through  a  theoretical  analysis 
of  the  towfish  dynamics  or  through  experiments,  then  a  greater  difference  may  result 
between  the  estimates  obtained  through  the  two  algorithms.  Notice,  however,  that  the 
confidence  intervals  resulting  from  the  extended  least-squares  algorithm  are  tighter  than 
those  obtained  through  deterministic  least-squares  estimation.  This  is  a  consequence 
of  the  reduction  of  the  variance  of  the  estimation  error  that  is  theoretically  expected 
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Figure  4.8:  Estimates  of  state-space  model  parameters  obtained  through  the  extended 
least-squares  algorithm  with  N  =  4  and  0  =  1.  These  figures  show  the  evolution  of  the 
estimates  of  the  AR  coefficients  for  Ax'(n]  (top),  Ay'[n]  (middle),  and  A0[n]  (bottom). 
For  simplicity,  the  cross-term  coefficients  are  not  shown. 
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to  result  from  the  incorporation  of  temporal  correlation  information  on  the  parameters 
through  the  state-space  model.  Figure  4.8  shows  the  estimates  of  the  AR  coefficients  for 
AxJ,[n],  Ay'[n]  and  A0[n].  For  simplicity,  the  cross-term  coefficients  (i.  e.,  those  relating 
Ai'[n]  to  A y'0[n)  and  A0[n],  etc.)  were  not  plotted. 

Figures  4.9  and  4.10  show  the  estimates  of  attitude  and  model  parameters  for  N  =  4 
and  0  =  0.987  and  the  same  initial  estimates  used  before.  In  this  case,  the  weight  of  the 
observations  decreases  progressively  as  they  become  “old,”  turning  the  algorithm  into  an 
adaptive  Kalman  filter  capable  of  responding  to  variations  in  the  dynamics  of  the  towfish 
and  in  the  statistics  of  the  observation  and  process  noises.  Taking  0  =  0.987  amounts 
to  applying  an  exponentially  decreasing  weight  to  the  data,  with  a  time  constant  of 
10  sec.  The  estimates  of  attitude  parameters  in  this  case  do  not  differ  much  from  those 
obtained  with  0  =  1.  However,  the  confidence  intervals  are  slightly  wider  because  the 
estimation  at  each  point  no  longer  takes  equally  into  account  all  previous  observations, 
but  gradually  “forgets”  old  observations.  The  estimates  of  the  AR  model  parameters  are 
shown  in  Fig.  4.10. 

The  estimates  of  Ai'[n],  Ay'[n]  and  A0[n]  obtained  in  this  chapter  will  be  used  in 
Chapter  5  for  correcting  the  geometric  distortions  in  the  sonograph. 
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Figure  4.10:  Estimates  of  state-space  model  parameters  obtained  through  the  extended 
least-squares  algorithm  with  N  =  4  and  ft  =  0.987.  These  figures  show  the  evolution 
of  the  estimates  of  the  AR  coefficients  for  Ai'[n]  (top),  A y'[n]  (middle),  and  A0[n] 
(bottom).  For  simplicity,  the  cross-term  coefficients  are  not  shown. 


96 


Chapter  5 


Correction  of  Geometric  Distortions 


The  distortion  parameters  estimated  through  the  technique  described  in  Chapter  4 
are  now  used  for  correcting  the  geometric  distortions  in  the  sonograph.  In  this  final  step 
the  sonograph  is  resampled  so  that  the  coordinates  of  a  point  in  the  processed  image 
correspond  more  closely  to  its  true  location  on  the  seabed,  assuming  the  estimates  of 
distortion  parameters  are  accurate.  Visual  inspection  of  the  corrected  sonograph  reveals 
considerable  reduction  of  the  geometric  distortions.  A  simulation  is  carried  out  to  provide 
further,  and  more  objective,  evidence  of  the  efficacy  of  our  technique  for  estimating  and 
correcting  geometric  distortions  in  sonographs. 


5.1  Correction  of  Cross-Track  Distortions 

As  indicated  in  Section  3.2,  the  slant-range  distortion  in  the  sonograph  has  to  be 
corrected  before  the  sampling  displacements  A nx'[/,  n]  can  be  estimated.  This  is  accom¬ 
plished  by  resampling  the  original  sonograph,  s[m,  n],  whose  m  coordinate  corresponds 
to  range  to  the  towfish,  to  obtain  a  new  image  s[/,  n],  whose  l  coordinate  corresponds  to 
horizontal  distances  on  the  seabed  plane. 

Let  s,(r,  n]  denote  the  sonograph  after  interpolation  along  each  line.1  The  slant-range 

lrrhe  combination  of  parenthesis  and  square  bracket  is  used  to  indicate  that  r  is  a  continuous  variable, 
while  n  is  a  discrete  variable. 
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Figure  5.1:  Height  of  the  water  column  measured  from  the  sonograph  after  lowpass 
filtering  to  eliminate  discontinuities  caused  by  quantization. 


corrected  image  is  given  by 


s[/,  n]  =  s,(r,  n 


r=v'P+(fc[nJ/(cT./2))a  ’ 


where  A[n]  is  the  distance  between  the  towfish  and  point  (ar0[nj,  y0[^]),  as  shown  in 
Fig.  2.4,  and  can  be  measured  from  the  sonograph  as  the  height  of  the  water  column 
at  each  line.  The  interpolation  is  carried  out  by  multiplying  the  discrete  Fourier  trans¬ 
form  of  the  line  by  a  phase  shift  corresponding  to  the  desired  sampling  location,  before 
calculating  the  inverse  transform.  However,  other  techniques,  such  as  truncated-sinc  or 
cubic-spline  interpolation,  may  be  used  instead.  If  S„[fc],  k  =  0, 1, ... ,  Nm  —  1  denotes 
the  DFT  of  the  nth  line  of  s[m,n],  we  have 

i[J,n]  =  S„[i]e*VP+W*<«r./a>)». 

k= 0 

Figure  5.1  shows  A[n]  for  the  sonograph  of  Fig.  1.2  after  lowpass  filtering  to  eliminate 
the  discontinuities  resulting  from  the  quantization  of  the  width  of  the  water  column  to 
an  integer  number  of  pixels.  The  slant-range  corrected  sonograph  was  shown  in  Fig.  3.4. 
From  it,  the  estimates  of  n]  are  calculated. 

As  discussed  in  Section  3.3,  before  the  correlation  lengths  can  be  calculated,  it  is 
necessary  to  shift  the  lines  of  the  slant-range  corrected  image,  s[/,n],  to  obtain  a  new 
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image,  ![/,  n],  for  which  the  sampling  displacements  A„x'[/,n]  are  zero.  According  to 
the  linear  model  of  Section  2.2.3,  at  each  line  of  the  image  we  have  A„x^[/,n]  =  ArJ,[n]. 
Therefore,  to  make  the  Anx'[/,  n]  equal  to  zero,  it  is  enough  to  shift  line  (n  +  1)  of  the 
image  by  —  Ax' [/,  n]/ (cT’,/2)  with  respect  to  line  n,  for  n  =  0, 1 , . . . ,  Arn  —  1 .  Equivalently, 
the  nth  line  of  the  image  has  to  be  shifted  by  Ax'0[i}/ (cTa/2)  with  respect  to  line 

number  zero.  That  can  be  done  by  multiplying  the  DFT  of  the  nth  line,  which  is  denoted 
by  Sn[fc],  by  the  appropriate  phase  shift,  and  then  taking  the  inverse  transform.  Thus, 

s[f,  n]  =  r-'  {£„[%-* ‘£"--  . 

The  estimates  Ax'[n]  used  in  this  operation  are  those  obtained  through  deterministic 
least-squares,  i.  e.,  they  are  simply  the  average  of  the  estimates  of  sampling  displacements, 
Anx' [/,  n]  for  each  line  of  the  image.  If  one  wishes  to  estimate  the  distortion  parameters 
through  adaptive  Kalman  filtering,  then  the  Ax' [n]  are  later  re-estimated  with  the  other 
two  distortion  parameters.  Ay' [n]  and  A0[n],  simultaneously,  through  the  extended  least- 
squares  algorithm. 

The  image  s[l,  n]  is  used  for  calculating  the  correlation  lengths,  L+[l,  n]  and  L~[l,n ], 
as  described  in  Section  3.3. 


5.2  Estimating  the  Sampling  Point  Coordinates 

After  the  distortion  parameters  are  estimated  through  one  of  the  techniques  pre¬ 
sented  in  Chapter  4,  the  slant-range  corrected  sonograph  is  resampled  for  correction  of 
the  remaining  geometric  distortions.  That  requires  obtaining  estimates  of  the  sampling 
coordinates,  x,[/,  n]  and  ya[l,  nj,  for  each  point  of  the  image,  from  the  estimates  of  dis¬ 
tortion  parameters,  Ax'[n],  Ay'[n]  and  A0[n].  Using  Eqs.  (2.12a)  and  (2.12b),  we  have 

A„x'[/,  n]  =  AxJ,[n]  (5.1a) 

AnJ/;[/,n]  =  Ay'0[n\  + l{cT,/2)A$[n\.  (5.1b) 

These  estimates  are  then  converted  from  the  (x',  y')  to  the  (x,  y)  coordinate  system.  From 
Eqs.  (2.13a)  and  (2.13b)  we  obtain 

x  =  (x'  —  x0[nj)  cos  0[n]  —  (y'  —  j/0[n])  sin  0[n] 
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y  =  (y'  -  y0[n])  cos  <?[n]  +  (x'  -  x0[n])  sin  0[r?]. 


Therefore  the  sampling  displacements  may  be  converted  from  one  coordinate  system  to 
the  other  through 

A  „x4[/,n]  =  A„x' [/,  n]  cos  0[n]  —  A„y'[/,  n]  sin  6[n\  (5.2a) 

A„y4[/,n]  =  Any'[/,  n]  cos  9[n]  +  A„ar' [/,  n]  sin  6[n\.  (5.2b) 

Finally,  the  estimates  of  sampling  point  coordinates  are  calculated  recursively  from  the 
estimates  of  sampling  displacements, 

x4[/,n  +  l]  =  x4[/,n]  +  A„x4[/,n]  (5.3a) 

y4[/,n+l]  =  y4[/,n]  +  A«y,[/,n],  (5.3b) 

starting  from  estimates  of  the  sampling  coordinates  of  the  first  line,  x4[/,  0]  and  y4[/,0], 
/  =  0,  ±1, . . . ,  zt(Ni  —  1).  Equations  (2.11a)  and  (2.11b)  may  be  used  to  calculate  these 
coordinates,  with  xo[0]  =  yo[0]  =  0  and  0[O]  =  <?[0]. 

A  graphic  representation  of  the  estimates  of  sampling  coordinates  can  be  provided  by 
plotting  the  estimated  sampling  points  (x4[m,  n],  y[m,  n])  in  a  mesh.  Figure  5.2  shows 
the  resulting  mesh  plot  for  our  sample  sonograph.  In  producing  this  figure  we  used  the 
estimates  of  distortion  parameters  obtained  through  deterministic  least-squares  estima¬ 
tion  in  Section  4.2  and  plotted  only  every  fourth  line  and  column  so  the  mesh  would 
not  be  too  dense.  This  plot  gives  a  clear  idea  of  how  the  sonar  beam  moved  acres*  the 
bottom  on  that  area,  according  to  the  estimates  of  distortion  parameters. 

5.3  Correcting  Geometric  Distortions 

We  now  address  the  issue  of  correcting  the  geometric  distortions  affecting  the  sono¬ 
graph.  The  procedure  is  divided  into  two  steps:  interpolation  in  the  cross-track  direction 
followed  by  interpolation  in  the  along-track  direction. 
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Figure  5.2:  Mesh  plot  illustrating  the  geometric  distortion  in  the  sonograph  of  Fig.  3.4. 
The  grid  marks  the  estimated  location  of  the  sampling  points  on  the  seabed  plane. 
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Figure  5.3:  Determining  the  points  for  interpolation  in  the  cross-track  direction.  Each 
line  of  the  image  is  resampled  at  points  that  are  equally  spaced  in  the  x  direction. 

5.3.1  Interpolation  in  the  Cross-TVack  Direction 

The  slant-range  corrected  image  s[/,  n]  is  first  interpolated  to  produce  a  new  image 
sc[l',n]  so  its  corresponding  sampling  points  on  the  seabed  are  equally  spaced  in  the  x 
direction.  Figure  5.3  illustrates  how  the  points  for  interpolating  each  line  of  s[/,n]  are 
determined.  The  black  dots  shown  in  the  figure  represent  the  estimates  of  the  location  of 
sampling  points  on  the  seabed  plane.  Because  of  the  linear  model  used  in  Chapter  2,  the 
estimates  of  sampling  points  for  each  sonar  pulse  fall  on  a  straight  line  that  approximates 
the  trajectory  of  the  sound  pulse  on  the  seabed  as  it  travels  away  from  the  towfish. 

The  sampling  points  of  sc[/',n]  are  indicated  in  the  image  by  the  uniformly  spaced 
circles  whose  cross-track  coordinates  are 

xc[V\  =  l\cTj 2),  /'  =  0,  ±1,  ±2, ... ,  ±{Nv  ~  1 ), 

where  Ni>  is  the  desired  number  of  columns  on  the  starboard  and  port  sides  of  the 
resampled  image.  The  corresponding  distances  along  the  nth  scan  line  are  denoted  by 
/c[/',  n]  and  in  general  are  non-integer.  They  are  calculated  by  finding  the  two  points 
[/x,n]  and  [/2, n]  whose  sampling  coordinate  estimates  x,[/j,n]  and  x,[/2,n]  are  located 
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immediately  before  and  immediately  after  £c[/'],  i.  e., 

i»[/i,n]  <  Xcll']  <  xt[l2,n],  l2  =  /i  +  1, 

and  then  by  computing  lc[l',  n]  by  linearly  interpolating  the  corresponding  /  coordinates, 

lc[l  ,n\  -  h  +  (*2  -  ‘1)7-77—1  .  f ,  V  (5-4) 

In  the  next  subsection,  we  will  need  to  know  the  along-track  coordinates  of  these  points, 
which  we  denote  by  yc[/',  n].  These  are  computed  by  linearly  interpolating  between  the 
y  coordinates  of  points  [/i,nj  and  [/2,n], 

STotC, r.]  =  -  ».['i."l)  HV[''17i‘[';;n|-r  (5.5) 

xa[lun\ 

The  nth  line  of  the  slant-range  corrected  image,  s[/,n],  is  then  interpolated  at  points 
/c[/',  n],  /'  =  0,1,...,  ( Ni>  —  1).  The  interpolation  is  again  performed  by  multiplying  the 
DFT  of  the  line  by  the  appropriate  phase  shift  for  each  sample  point,  before  calculating 
the  inverse  transform.  Thus,  if  5„[fc],  k  =  0, 1, . . . ,  JVj  —  1  denotes  the  DFT  of  the  n-th 
line  of  s[/,n],  we  have 

Jt=o 

5.3.2  Interpolation  in  the  Along-Track  Direction 

By  following  the  procedure  of  the  last  subsection  one  obtains  a  new  image  sc[l',  n], 
whose  corresponding  sampling  points  on  the  seabed  are  now  correctly  located  in  the  x 
direction,  but  are  still  irregularly  spaced  in  the  y  direction.  The  final  step  for  correcting 
the  geometric  distortions  is  to  resample  each  column  of  sc[l',  n]  to  produce  a  final  image 
sa[l\  n']  whose  corresponding  sampling  points  on  the  seabed  plane  are  regularly  spaced 
in  both  the  x  and  y  directions. 

Figure  5.4  illustrates  this  new  situation.  The  circles  shown  in  that  figure  are  the  same 
as  those  shown  in  Fig.  5.3  and  mark  the  location  of  the  sampling  points  associated  with 
sc[V,  n];  the  points  where  each  column  should  be  interpolated  are  marked  in  the  figure 
by  small  boxes.  Their  along-track  coordinates  are  given  by 

t/aM  =  n'(v7»,  n'  =  0,1,..., 7V„.  -  1, 
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Figure  5.4:  Determining  the  points  for  interpolation  in  the  along-track  direction. 

where  jVn<  is  the  desired  number  of  lines  in  the  final  image.  The  corresponding  distances 
along  the  /'- th  column  of  sc[/',  n]  are  denoted  by  na[l',n'}  and  are  in  general  non-integer. 
Paralleling  the  procedure  of  the  previous  section,  these  are  calculated  by  finding  the  two 
points  [/', nj]  and  [/', n2]  whose  sampling  coordinates  yc[l',n  1]  and  yc[/', «2],  calculated 
through  Eq.  (5.5)  are  located  immediately  before  and  immediately  after  ya[n'],  i.  e., 

1]  <  S/a[«']  <  yc[l\n2},  n2  =  nx  +  1, 


and  then  by  computing  na[l',  n'j  by  linearly  interpolating  the  corresponding  n  coordinates, 


=  nj  +  ( n2  -  nj) 


~  yc[l',ni] 

yc[l\n2]  -  yc[l',Tii] 


Before  proceeding  to  interpolate  the  columns,  we  need  to  address  a  problem  that  is 
not  present  in  the  cross-track  direction.  Because  of  the  occurrence  of  backscanning,  the 
sequence  of  coordinates  t/c[/',  n],  n  =  0, 1, . . . ,  Nn  —  1,  may  not  be  strictly  increasing,  i.  e., 
we  do  not  necessarily  have  yc[/',  ni]  >  yc[l',n 2]  for  ri\  >  n2.  As  pointed  out  in  Chapter  1, 
backscanned  areas  appear  in  triplicate  in  the  image,  although  the  three  passes  may 
not  look  exactly  the  same  due  to  the  different  viewing  angles.  Because  of  that  fact,  and 
because  the  registration  among  the  three  passes  may  not  be  accurate  enough,  attempting 


104 


to  merge  the  three  views  of  a  backscanned  area  during  the  interpolation  of  the  columns 
is  not  likely  to  yield  good  results.  Instead,  we  choose  to  preserve  only  one  view  of  the 
backscanned  area.  To  achieve  greater  image  continuity,  it  is  better  to  select  either  the 
first  or  the  third  pass,  since  the  middle  pass  is  isolated  from  both  the  preceding  and 
the  forthcoming  portions  of  the  sonograph.  Furthermore,  the  beam  is  likely  to  scan  the 
bottom  more  slowly  during  the  first  pass,  when  it  is  de-accelerating  to  reverse  its  direction 
of  movement,  than  during  the  third  pass  when,  after  reversing  direction  a  second  time,  it 
is  accelerating  to  compensate  for  the  time  lost  during  backscanning.  Therefore,  we  select 
only  the  first  pass  over  backscanned  areas.  In  practice,  that  amounts  to  considering  only 
progressively  increasing  values  of  yc[l',  n]  when  looking  for  the  interval  surrounding  the 
interpolating  coordinate  ya[n']. 

Resampling  in  the  along-track  direction  is  again  accomplished  by  multiplying  the 
DFT  of  each  column  of  sc[l,  n']  by  the  appropriate  phase  shift  for  each  point  before 
computing  the  inverse  transform.  Denoting  by  k  =  0,1,...,  Nn  —  1,  the  DFT  of 

the  /' th  column  of  sc[/',  n],  we  have 

k= 0 

Figure  5.5  shows  the  sonograph  of  Fig.  1.2  after  correction  of  the  geometric  distortions 
as  described  in  this  section.  The  estimates  of  distortion  parameters  used  in  estimating  the 
sampling  coordinates  are  those  obtained  through  deterministic  least-squares  estimation 
in  Chapter  4.  The  lower  right  corner  of  the  image  is  shown  enlarged  in  Fig.  5.6  alongside 
the  corresponding  area  of  the  original  sonograph.  Notice  how  the  cable  lying  on  the 
bottom  that  appeared  pronouncedly  jagged  as  a  result  of  yawing  now  presents  a  slowly 
varying  curvature  that  is  likely  closer  to  its  true  shape.  Notice  also  how  the  multiple 
images  of  objects  in  backscanned  areas  are  correctly  replaced  by  single  images  of  those 
objects.  However,  some  areas  of  the  corrected  image  in  Fig.  5.5  appear  slightly  blurred, 
as  a  result  of  aliasing  during  the  acquisition  of  the  original  sonograph.  In  fact,  perfect 
reconstruction  of  the  image  is  not  possible  in  areas  where  the  sampling  rate  in  the  along- 
track  direction  is  greater  than  the  Nyquist  rate  of  the  backscattering  function  b(x,y).  as 
a  result  of  faster  motion  of  the  sonar  beam. 
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Figure  5.5:  Sonograph  of  Fig.  1.2  after  correction  of  geometric  distortions. 


Figure  5.6:  Comparison  of  a  detail  of  the  original  and  corrected  sonographs.  Correction 
of  geometric  distortions  eliminates  the  jagged  appearance  of  the  cable  and  the  multiple 
images  of  rocks,  caused  by  barkscanning. 
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5.4  Simulation 


The  side-scan  sonar  employed  for  making  the  sonographs  used  in  this  thesis  was  a  common 
commercial  unit  and  was  not  equipped  with  attitude  sensors.  Therefore,  measurements 
of  the  true  values  of  the  attitude  parameters  are  not  available  for  comparison  with  the 
estimates  obtained  in  Chapter  4.  As  seen  in  Section5.3,  visual  inspection  of  the  sonograph 
after  correction  of  the  geometric  distortions  provides  strong  subjective  indication  of  the 
efficacy  of  our  algorithm.  However,  a  more  objective  evaluation  can  be  provided  by  a 
simulation. 

Simulating  geometric  distortions  requires  an  undistorted  image  of  the  seabed  that 
can  be  resampled  in  an  irregular  pattern  to  simulate  the  effect  of  towfish  instabilities. 
Since  all  the  sonographs  in  our  data  set  present  some  degree  of  geometric  distortion,  it 
is  necessary  to  artificially  generate  an  image  on  which  to  perform  the  simulation.  The 
sonograph  of  Fig.  5.5  itself  could  be  used  to  that  end,  provided  that  application  of  the 
algorithm  of  Chapters  3  and  4  to  that  corrected  image  indicates  that  it  is  free  of  geometric 
distortions.  However,  such  is  not  the  case,  because  of  the  occurrence  of  aliasing  in  the 
original  sonograph.  As  explained  in  Section  5.3,  in  areas  where  the  original  sonograph 
is  aliased,  the  corrected  image  appears  blurred,  which  causes  the  correlation  lengths  to 
be  longer  than  in  areas  free  of  aliasing.  As  a  result,  the  algorithm  indicates  that  the 
corrected  image  is  still  affected  by  geometric  distortions. 

To  obtain  the  desired  image  for  the  simulation,  it  is  necessary  to  circumvent  the 
problem  presented  by  aliased  areas.  In  Chapter  4  the  distortion  parameters  for  each 
line  of  the  sonograph  were  estimated  from  the  side  that  had  smaller  estimated  sampling 
intervals  A„y'[/,  n].  That  situation  was  illustrated  in  Fig.  4.1.  In  that  example,  when  the 
geometric  distortions  are  corrected  the  sampling  displacements  on  the  port  side  (/  <  0) 
are  made  larger  to  fit  the  linear  model  represented  by  the  straight  line.  That  results  in 
the  corrected  image  appearing  blurred  on  the  port  side  at  that  line,  as  explained  above. 
That  effect  can  be  avoided  if  the  distortion  parameters  are  estimated  separately  from 
each  side  of  the  line,  as  indicated  in  Fig.  5.7.  Figure  5.8  shows  the  result  of  processing 
each  side  of  the  sonograph  with  estimates  of  distortion  parameters  obtained  from  that 
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Figure  5.7:  An  image  unaffected  by  aliasing-related  blurring  can  be  obtained  by  estimat¬ 
ing  the  distortion  parameters  separately  from  each  side  of  the  sonograph,  as  illustrated 
in  this  example. 

side  only.  The  rectangular  areas  outlined  in  each  image  are  then  cut  out  and  spliced  side 
by  side  to  produce  the  image  in  Fig.  5.9,  that  can  be  used  for  the  simulation. 

The  image  of  Fig.  5.9  is  now  resampled  on  an  irregular  grid  to  simulate  the  production 
of  a  sonograph  affected  by  geometric  distortions.  The  simulated  distortion  parameters 
are  generated  by  passing  white  Gaussian  noise  through  an  autoregressive  model  with 
4  poles  of  magnitude  0.8  and  phases  ±?r/8  and  ±ir/16  radians.  The  location  of  the 
poles  and  the  energy  of  the  input  noise  were  chosen  so  that  the  resulting  simulated 
parameters  had  characteristics  similar  to  those  of  the  distortion  parameters  estimated 
from  the  original  sonograph.  Figure  5.10  shows  the  undistorted  image  and  the  simulated 
distorted  sonograph.  The  figure  also  shows  the  simulated  sonograph  after  correction  of 
the  slant-range  distortion  and  detection  of  backscanning. 

Figure  5.11  shows  the  simulated  attitude  parameters,  along  with  the  estimates  ob¬ 
tained  through  deterministic  least-squares  estimation,  adaptive  Kalman  filtering,  and 
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Figure  5.8:  Result  of  independently  processing  the  starboard  and  port  sides  of  the 
graph,  to  avoid  blurring  of  aliased  areas. 


Figure  5.9:  Image  used  in  the  simulation.  It  was  formed  by  joining  the  rectangular 
sections  outlined  in  the  images  of  Fig.  5.8. 

with  a  Kalman  filter  using  the  true  state-space  model  employed  to  generate  the  simu¬ 
lated  distortion  parameters.  The  estimates  obtained  through  adaptive  Kalman  filtering 
do  not  differ  significantly  from  those  obtained  through  deterministic  least-squares  estima¬ 
tion,  except  for  a  smaller  high-frequency  content.  As  noted  in  Chapter  4,  this  is  because 
all  the  autoregressive  parameters  are  being  estimated  concurrently  with  the  distortion 
parameters  and,  consequently,  the  state-space  model  contributes  little  a  priori  informa¬ 
tion.  In  fact,  the  estimates  of  state-space  model  parameters,  a[n],  obtained  through 
the  adaptive  Kalman  filter  failed  to  converge  to  the  true  values.  However,  if  the  true 
state-space  model  is  used  in  a  standard  Kalman  filter,  the  resulting  estimates  show  a 
noticeable  improvement,  as  seen  in  the  figure. 

Table  5.1  presents  the  maximum  absolute  estimation  errors  for  the  deterministic 
least-squares  (DLS),  adaptive  Kalman  filter  (AKF),  and  standard  Kalman  filter  (I\F) 
algorithms.  Notice  that  using  the  true  state-space  model  in  the  standard  Kalman  filter 
reduces  the  maximum  absolute  estimation  error  of  all  three  parameters  by  approximately 
50  %  with  respect  to  deterministic  least-squares  estimation  and  adaptive  Kalman  filtering. 
The  estimates  of  yaw  angle  obtained  through  any  of  the  three  algorithms  are  remarkably 
accurate,  with  the  standard  Kalman  filter  yielding  a  maximum  absolute  estimation  er¬ 
ror  of  less  than  one  degree.  The  estimates  presented  here  were  obtained  with  p  =  140. 
Smaller  estimation  errors  may  be  expected  from  all  algorithms  with  a  larger  number  of 
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Figure  5.10:  The  area  outlined  in  the  undistorted  image  (top)  is  iesampled  as  if  scanned 
by  a  side-scan  sonar  subject  to  motion  instabilities,  to  produce  a  simulated  distorted 
sonograph  (middle).  The  bottom  image  is  the  simulated  sonograph  after  correction  of 
the  slant-range  distortion  and  detection  of  backscanning. 
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Figure  5.11:  Simulated  attitude  parameters  (SIM)  and  estimates  obtained  through  deter¬ 
ministic  least-squares  estimation  (DLS)  ,  adaptive  Kalman  filtering  (AKF),  and  standard 
Kalman  filtering  (KF). 
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Table  5.1:  Maximum  absolute  estimation  errors. 


Algorithm 

Xf 

<t> 

e 

DLS 

26.87  cm 

11.37° 

1.18° 

AKF 

29.25  cm 

11.03° 

1.19° 

KF 

18.84  cm 

4.84° 

0.66° 

observation  points. 

The  results  obtained  with  the  adaptive  Kalman  filter  and  with  the  standard  Kalman 
filter  correspond  to  extreme  cases  in  which  either  none  of  the  autoregressive  model  param¬ 
eters  are  known  or  all  of  them  are  known,  respectively.  Intermediate  values  of  maximum 
estimation  errors,  as  well  as  better  estimates  of  the  unknown  model  parameters,  can 
be  expected  from  the  adaptive  Kalman  filter  if  only  a  few  of  the  autoregressive  model 
parameters  are  unknown. 

Figure  5.12  shows  the  undistorted  image  along  with  the  simulated  sonograph  and  the 
reconstructed  image.  Comparison  of  the  shapes  of  rocks  and  scours  in  the  three  images 
shows  that  the  algorithm  achieves  a  significant  reduction  of  the  simulated  geometric 
distortions.  As  pointed  out  before,  perfect  reconstruction  is  not  possible  in  areas  where 
the  sampling  interval  was  large  enough  to  cause  aliasing. 


5.5  Further  Examples 

A  second  sonograph  from  our  data  set  is  shown  in  Fig.  5.13,  and  in  Fig.  5.14  after  correc¬ 
tion  of  the  slant-range  distortion  and  backscanning  detection.  The  corrected  sonograph 
is  shown  in  Fig.  5.15.  The  correction  of  geometric  distortions  is  most  evident  in  the  up¬ 
per  right-hand  corner  of  the  image,  where  features  on  the  bottom  appear  pronouncedly 
stretched  in  the  original  sonograph,  but  acquire  more  natural  shapes  after  processing. 
The  estimates  of  attitude  parameters  obtained  from  that  sonograph  through  determinis¬ 
tic  least-squares  estimation  are  shown  in  Fig.  5.16. 

A  third  example  is  given  in  Figs.  5.17  through  5.19.  As  in  the  first  two  examples, 
the  jagged  appearance  of  the  cable  seen  in  the  original  sonograph  is  made  consider- 
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Figure  5.12:  Comparison  of  the  original  image  (topj.  t he  simulated  sonograph  (middle) 
and  reconstructed  image  (bottom)  reveals  that  the  algorithm  is  capable  of  effectively 
reducing  geometric  distortions. 
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Figure  5.13:  A  second  sonograph  from  our  data  set. 
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backscanning  detection. 


Figure  5.15:  Sonograph  of  Fig.  5.13  after  correction  of  geometric  distortions. 


Figure  5.16:  Estimates  of  attitude  parameters  for  the  sonograph  of  Fig.  5.13. 
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ably  smoother  by  the  correction  of  geometric  distortions.  Several  rocks  that  appear 
compressed,  stretched,  or  with  multiple  images  in  the  original  sonograph  acquire  more 
natural  shapes  after  processing,  which  are  likely  closer  to  their  true  shapes  on  the  seabed. 

Geometrically  corrected  sonographs,  such  as  those  shown  in  Figs.  5.5,  5.15,  and  5.19, 
are  the  end  product  of  the  techniques  presented  in  this  thesis.  The  next  chapter  presents 
final  considerations. 
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Figure  5.17:  A  third  sonograph  from  our  data  set. 


Figure  5.18:  Sonograph  of  Fig.  5.17  after  correction  of  the  slant-range  distortion  and 
backscanning  detection. 
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Figure  5.19:  Sonograph  of  Fig.  5.17  after  correction  of  geometric  distortions. 
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Figure  5.20:  Estimates  of  attitude  parameters  for  the  sonograph  of  Fig.  5.17. 
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Chapter  6 


Final  Considerations 


In  this  chapter  we  offer  final  considerations  concerning  the  theoretical  developments, 
algorithms  and  results  presented  in  this  thesis,  along  with  suggestions  for  future  research 
on  the  correction  of  geometric  distortions  in  sonographs. 


6.1  Evaluation  of  Results 

Visual  inspection  of  the  corrected  sonographs  presented  in  Chapter  5  offers  strong 
indication  that  the  techniques  developed  in  this  thesis  can  effectively  reduce  the  degree 
of  geometric  distortion  in  sonographs,  thus  fulfilling  their  main  goal.  The  algorithm  is 
original  in  that  it  does  not  require  navigational  or  attitude  measurements  for  correcting 
the  geometric  distortions,  but  relies  solely  on  the  information  contained  in  the  image  itself, 
and  on  certain  assumptions  about  its  statistical  properties.  In  the  process  of  correcting 
the  geometric  distortions,  the  algorithm  produces  other  relevant  information  that  can 
prove  useful  in  the  interpretation  of  sonographs.  For  instance,  it  provides  a  means  of 
detecting  lines  affected  by  backscanning,  apparently  with  relatively  low  percentages  of 
false  alarms  and  misses,  as  seen  in  Chapter  3.  It  also  allows  reconstruction  of  the  sampling 
pattern  on  the  bottom  in  the  form  of  mesh  plots  such  as  the  one  presented  in  Chapter  5. 

Another  by-product  of  the  algorithm  are  estimates  of  the  towfish  attitude  parame¬ 
ters,  namely  its  lateral  displacement  in  the  cross-track  direction  and  its  pitch  and  yaw 
angles.  Because  of  the  accumulation  of  the  estimation  errors  of  parameter  increments, 
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the  estimates  of  attitude  parameters  become  progressively  less  reliable,  as  indicated  by 
the  widening  of  the  confidence  intervals  shown  in  the  plots  of  Chapter  4,  and,  conse¬ 
quently,  so  do  the  estimates  of  sampling  point  coordinates.  The  estimates  of  distortion 
parameters,  Ax'0 [n],  Ay'0[n}  and  A0[n],  however,  remain  just  as  accurate,  which  means 
that  the  relative  location  of  points  of  the  image  with  respect  to  neighboring  points  can 
be  estimated  more  accurately  than  their  absolute  location.  In  other  words,  the  algorithm 
is  more  capable  of  correcting  local  geometric  distortion  than  of  determining  the  absolute 
location  of  objects  in  the  image. 

The  simulation  carried  out  in  Chapter  5  provides  more  objective  evidence  of  the  effi¬ 
cacy  of  the  different  forms  of  the  algorithm.  According  to  that  simulation,  very  accurate 
estimates  of  the  towfish  yaw  angle  are  obtained  through  both  deterministic  least-squares 
estimation  and  adaptive  Kalman  filtering,  as  well  as  fairly  accurate  estimates  of  its  hor¬ 
izontal  lateral  displacement.  The  estimates  of  pitch  are  less  reliable,  because  variations 
in  pitch  angle  produce  less  geometric  distortion  than  equal  variations  in  yaw  angle. 

The  estimates  of  distortion  and  attitude  parameters  obtained  through  the  adaptive 
Kalman  filter  do  not  differ  considerably  from  those  obtained  through  deterministic  least- 
squares  estimation,  because  all  the  parameters  of  the  state-space  model  used  in  the  filter 
are  unknown  and  are  simultaneously  estimated  with  the  distortion  parameters.  In  that 
case,  the  state-space  model  contributes  little  a  priori  information.  Improved  results  can 
be  expected  if  the  state-space  model  is  derived  from  a  theoretical  analysis  of  the  dynamics 
of  the  towfish  and  only  a  few  of  its  parameters  are  unknown. 

6.2  Suggestions  for  Future  Research 

The  most  important  path  for  future  research  on  the  technique  presented  in  this 
thesis  is  a  more  definitive  evaluation  of  its  accuracy.  One  way  that  may  be  accomplished 
is  through  further  and  more  extensive  simulations.  Definitive  evaluation,  however,  can  be 
accomplished  by  comparing  the  estimates  of  attitude  parameters  with  real  measurements 
obtained  through  attitude  sensors  mounted  on  a  towfish.  Another  possible  means  of 
evaluation  is  to  acquire  sonographs  of  areas  for  which  sufficiently  detailed  bathymetric 
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charts  are  available,  and  then  attempt  to  match  points  in  the  corrected  sonographs  with 
points  in  the  chart. 

Another  goal  for  future  research  could  be  to  make  the  estimates  of  distortion  parame¬ 
ters  more  robust  by  improving  the  process  of  selection  of  observation  points  for  each  line 
of  the  sonograph.  As  mentioned  in  Chapter  3,  that  would  probably  require  the  inclusion 
of  image  segmentation  techniques  to  divide  the  sonograph  in  areas  where  the  correlation 
lengths  may  be  assumed  constant.  A  more  straightforward  modification  that  can  also  be 
made  to  the  algorithm  is  to  incorporate  a  high-pass  filter  or  some  other  device  to  period¬ 
ically  correct  the  estimates  of  attitude  parameters  in  order  to  keep  them  from  diverging 
as  a  result  of  the  accumulation  of  the  estimation  errors  of  the  parameter  increments. 

As  pointed  out  in  Chapter  4,  our  intent  in  employing  a  Kalman  filter  in  the  estimation 
of  distortion  parameters  is  to  provide  a  framework  in  which  knowledge  of  the  dynamics  of 
the  towfish  may  be  readily  incorporated  in  the  future.  Improvements  in  the  estimates  can 
in  principle  be  obtained  by  replacing  the  general  auto-regressive  model  we  adopted  with 
a  more  specific  state-space  model  derived  from  a  mechanical  analysis  of  the  dynamics  of 
towed  cylinders.  Additional  improvements  can  be  expected  if  the  estimation  of  distortion 
parameters  is  carried  out  through  smoothing  or  fixed-lag  smoothing  instead  of  through 
regular  Kalman  filtering,  and  if  a  greater  number  of  observation  points  is  used. 

On  a  more  practical  side,  an  issue  to  be  investigated  is  the  feasibility  of  implementing 
the  algorithm  in  real  time,  so  that  the  corrected  image  may  be  seen  on  a  display  in  the 
deploying  vessel  as  the  towfish  scans  the  bottom. 


6.3  Conclusions 

The  main  contribution  of  this  thesis  is  the  development  of  a  technique  for  the  estima¬ 
tion  and  correction  of  geometric  distortions  in  side-scan  sonar  images.  Other  techniques 
previously  reported  in  the  literature  utilize  navigational  data  to  correct  large-scale  dis¬ 
tortions  caused  by  variations  in  the  course  of  the  deploying  vessel  and  by  the  slant-range 
efFect.  In  the  case  of  a  relatively  few  more  sophisticated  units  equipped  with  attitude  sen¬ 
sors,  it  is  also  possible  to  correct  geometric  distortions  due  to  motion  instabilities  of  the 


127 


towfish.  By  comparison,  in  the  technique  presented  in  this  thesis  the  geometric  distortion 
is  estimated  from  the  image  itself,  requiring  no  navigational  or  attitude  measurements.  It 
may  thus  be  applied  to  any  existing  digitized  sonographs  or  to  analog  sonographs  stored 
on  magnetic  tape. 

A  very  detailed  analysis  of  the  side-scan  sonar  geometry  was  carried  out,  leading  to  a 
nonlinear  model  that  expresses  the  sampling  point  coordinates  as  a  function  of  the  towfish 
attitude  parameters.  Unlike  previous  studies  found  in  the  side-scan  sonar  literature  (see 
Chapter  2),  this  model  takes  into  account  the  misalignment  of  the  transmitting  and 
receiving  beams. 

The  techniques  presented  here  are  likely  to  find  wider  application  in  the  enhancement 
of  high-resolution  sonographs  (roughly,  those  with  operating  frequencies  of  100  kHz  or 
higher).  High-resolution  side-scan  sonars  are  usually  very  compact,  portable  units  in 
which  the  transducer  array  is  mounted  on  a  relatively  small  towfish.  In  this  type  of  sonar, 
the  towfish  is  typically  not  equipped  with  sensors  to  provide  attitude  measurements  that 
could  be  used  for  correcting  geometric  distortions.  Furthermore,  their  small  size  makes 
them  more  susceptible  to  motion  instabilities  than  their  larger  and  heavier  long-range 
counterparts.  Those  factors  make  this  type  of  side-scan  sonar  the  most  likely  beneficiary 
of  the  techniques  presented  here. 
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Appendix  A 


Summary  of  Notation 


The  tables  presented  in  this  Appendix  list  the  definitions  of  variables  by  chapter,  in 
order  of  appearance.  For  the  complete  definition  of  each  variable,  refer  to  the  sections, 
equations  or  figures  listed  in  the  reference  columns. 


Table  A.l:  Definitions  from  Chapter  2 


Symbol 

Description 

* 

Backscattering  strength  on  the  seabed. 

Sec.  2.1 

s 

Image  intensity  of  the  digitized  sonograph. 

Sec.  2.1 

[m,r?] 

Coordinate  system  defined  on  the  sonograph. 

Fig.2.1 

Coordinate  system  defined  on  the  seabed. 

Fig-  2.1 

ixsi  y») 

Coordinates  of  sampling  points. 

Sec.  2.1,  Fig.  2.4 

(*j,yj,z/) 

Coordinates  of  the  towfish. 

Fig.  2.3 

4> 

Pitch  angle  of  the  towfish. 

Figs.  1.4,  2.3 

e 

Yaw  angle  of  the  towfish. 

Figs.  1.4,  2.3 

t 

Time  from  start  of  acquisition  of  the  sonograph. 

Sec.  2.2.1 

tj 

Time  when  a  pulse  is  transmitted. 

Sec.  2.2.1 

tR 

Time  when  the  returned  signal  is  sampled. 

Sec.  2.2.1 

{xe,ye,ze) 

Point  on  the  axial  plane  of  the  effective  beam. 

Eq.  (2.3),  Fig.  2.4 

4>e,Qe 

Pitch  and  yaw  angle  of  the  effective  beam. 

Eq.  (2.2),  Fig.  2.4 

( xo,y0 ) 

Reference  point  in  the  effective  beam  plane. 

Fig.  2.4 

d 

Distance  between  ( x$,yg )  and  (x0,y0). 

Fig.2.3 

h 

Distance  between  {x0,y0)  and  the  towfish. 

Fig.  2.3 

r 

Range  from  a  sampling  point  to  the  towfish. 

Sec.  2.2.1 

c 

Speed  of  sound  in  water. 

Sec.  2.2.1 
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Table  A. 2:  Definitions  from  Chapter  2  (Continued) 


Symbol 

Description 

Reference 

Tj 

T , 

Am 

Nn 

S 

N, 

AnXa 

A„j/a 

(*',*') 

Firing  period  of  the  sonar. 

Sampling  period  used  in  digitizing  the  signal. 
Number  of  samples  per  line  of  the  sonograph. 
Number  of  lines  in  the  sonograph. 

Image  after  correction  of  the  slant-range  distortion. 
Number  of  pixels  per  line  per  channel  of  s. 

Sampling  displacement  in  the  x  direction. 

Sampling  displacement  in  the  y  direction. 

Auxiliary  coordinate  system. 

Sec.  2.2.1 

Sec.  2.2.1 

Sec.  2.2.1 

Sec.  2.2.1 

Sec.  2.2.2 

Sec.  2.2.2 

Eq.  (2.7a),  Fig.  2.5 
Eq.  (2.7b),  Fig.  2.5 
Fig.  2.6 

Table  A. 3:  Definitions  from  Chapter  3 


Symbol 

Description 

Reference 

Rb 

Autocorrelation  function  of  b(x,y). 

Sec.  3.1 

Upsampled  sonograph  image. 

Sec.  3.1 

Bf  '"dt 

Upsampling  factor. 

Sec.  3.1 

HP 

Cross-track  coordinate  in  the  upsampled  image. 

Sec.  3.1 

I'.". 

Along-track  coordinate  in  the  upsampled  image. 

Sec.  3.1 

Mean  of  line  segments. 

Eq.  (3.1) 

wm 

Standard  deviation  of  line  segments. 

Eq.  (3.2) 

p 

Normalized  cross-correlation  of  line  segments. 

Eq.  3.3,  Fig.  3.1 

A  V 

Relative  shift  between  line  segments  in  the  V  direction. 

Fig.  3.1 

An' 

Relative  shift  between  line  segments  in  the  n'  direction. 

Fig.  3.1 

21  +  1 

Length  of  line  segments  for  calculating  cross-correlations. 

Fig.  3.1 

Pn 

Sampling  interval  in  the  along-track  direction. 

Fig.  3.9 

Positive-  and  negative-lag  correlation  lengths. 

Fig.  3.7 

Po 

Threshold  for  calculating  the  correlation  length. 

Fig.  3.7 

U 

Correlation  length  of  the  backscattering  function  6(x,y). 

Sec.  3.3.1 

Ljio 

Average  correlation  length  in  the  absence  of  distortions. 

Sec.  3.3.1 

9 

Scaled  sampling  displacements. 

Eq.  (3.9) 

D 

Matrix  of  coefficients  for  the  system  of  equations  on  q. 

Sec.  3.3.1 

h 

Inhomogeneous  vector  for  the  system  of  equations  on  q. 

Eq.  (3.11) 

X 

Set  of  lines  chosen  for  estimation  of  LTO. 

Sec.  3.8 

Nm 

Number  of  lines  in  M. 

Sec.  3.8 

V 

Average  speed  of  the  deploying  vessel. 

Sec.  3.8 

Average  speed  of  the  towfish  over  the  lines  in  N. 

Sec.  3.8 

130 


Table  A. 4:  Definitions  from  Chapter  4 


Symbol 

Description 

Reference 

C 

Set  of  distances  of  selected  observation  points. 

Sec.  4.1 

P 

Number  of  points  in  C 

Sec.  4.1 

cx?  cy 

Chi-square  cost  functions. 

Sec.  4.2 

•b 

Variance  of  the  estimation  error  of  Ai' 

Sec.  4.2 

atl 

At/1 

Variance  of  the  estimation  error  of  A y'0 

Sec.  4.2 

Variance  of  the  estimation  error  of  A 6 

Sec.  4.2 

4. 

Variance  of  the  estimation  error  of  xa 

Sec.  4.2 

< 

Variance  of  the  estimation  error  of  y0 

Sec.  4.2 

Variance  of  the  estimation  error  of  6 

Sec.  4.2 

c 

Observation  vector. 

Eq.  (4.13) 

i 

Vector  of  distortion  parameters. 

Eq.  (4.13) 

£ 

Observation  noise  vector. 

Eq.  (4.13) 

c 

Observation  matrix. 

Eq.  (4.14) 

A<, 

Covariance  matrices  of  £  and 

Sec.  4.3.1 

Cross-covariance  matrices  of  £  with  C  and  v. 

Sec.  4.3.1 

N 

Order  of  the  state-space  model. 

Sec.  4.3.2 

kmL 

State  vector. 

Sec.  4.3.2 

A 

State  transition  matrix. 

Eq.  (4.16) 

B 

Input  matrix. 

Eq.  (4.16) 

C 

Redefined  observation  matrix. 

Eq.  (4.16) 

in 

Process  noise  vector. 

Eq.  (4.16) 

R 

Covariance  matrix  of  ui- 

Sec.  4.3.2 

S 

Cross-covariance  matrix  of  v  and  tv. 

Sec.  4.3.2 

& 

Auxiliary  vector  of  the  information  filter. 

Sec.  4.3.3 

P 

Covariance  matrix  of  the  state  estimation  error. 

Sec.  4.3.3 

G,N 

Auxiliary  matrices  in  the  information  filter. 

Sec.  4.3.3 

0 

Forgetting  factor. 

Secs.  4.3.3,  4.3.4 

Ai,...,Ah 

Matrices  of  the  multi-dimensional  AR  model. 

Sec.  4.3.4 

a 

Vector  of  model  parameters. 

Sec.  4.3.4 

X 

Regression  matrix. 

Sec.  4.3.4 

U 

Covariance  matrix  of  the  parameter  estimation  error. 

Sec.  4.3.4 

J 

Auxiliary  matrix  in  the  RLS  algorithm. 

Sec.  4.3.4 
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Table  A. 5:  Definitions  from  Chapter  5 


Symbol 

Description 

Reference 

sc 

Image  after  correction  in  the  cross-track  direction. 

Sec.  5.2.1 

V 

Cross-track  coordinate  of  $c. 

Sec.  5.2.1 

Nr 

Number  of  pixels  per  line  per  channel  in  sc. 

Sec.  5.2.1 

lc 

Coordinates  for  resampling  the  lines  of  s. 

Sec.  5.2.1 

(*c,ye) 

Coordinates  of  sampling  points  associated  with  sc. 

Sec.  5.2.1 

Sn 

DFT  of  the  n-th  line  of  s 

Sec.  5.2.1 

sa 

Final  image  after  correction  of  geometric  distortions. 

Sec.  5.2.2 

n' 

A  long-track  coordinate  of  sa. 

Sec.  5.2.2 

Nn> 

Number  of  pixels  per  column  in  sa. 

Sec.  5.2.2 

na 

Coordinates  for  resampling  the  columns  of  sc. 

Sec.  5.2.2 

Sr 

DFT  of  the  /'-th  line  of  sc 

Sec.  5.2.2 
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Appendix  B 


Derivation  of  Results  from  the 
Geometry  of  the  Problem 


In  Section  2.2.1  it  was  necessary  to  express  the  distance  d  defined  in  Fig.  2.4  as  a  function 
of  the  towfish  attitude  parameters  at  times  tj  and  tp.  The  key  to  deriving  this  expression 
is  that  the  round-trip  travel  time  of  the  wavefront  from  the  point  where  the  towfish  was 
located  at  time  tj  to  the  sampling  point  (xt(tp),y,(tp))  and  back  to  the  point  where 
the  towfish  was  located  at  time  tp  has  to  equal  ( tp  —  tr).  Let  us  denote  by  rj  and 
rp  the  distances  from  points  (xj{tT),yj(tT),zj{tT))  and  (xj(tp),  yf(tp),  z/(tp))  to  point 
(x,{tp),yt(tR)),  i-  e., 

rT  =  ([x/(*r)  ~  +  [y/(*r)  -  !/#(<k)]2  +  z/(*r))2  (B.l) 

rR  =  ([*/(*«)  -  +  [y/(<*)  -  y*(**)]2  +  z/(t r )) 2  •  (B.2) 

Then,  as  argued  above,  we  must  have 


rT  +  rR  =  c(tp  —  tr), 


where  c  is  the  speed  of  sound  in  water.  This  equation  describes  a  three-dimensional 
ellipsoid  with  focii  at  points  (x/(tr),  j//(<r)i  */(<r))  and  (x/(tp),y/(tp),  Zf(tp)).  We  will 
denote  by  r  the  average  of  rj  and  rp, 


A  rj  +  rp 

r  = - 

2 


(B.3) 
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Thus,  r  is  also  given  by 


c(tn  —  t  t) 
r~  2 

From  Fig.  2.4  we  had  obtained  Eqs.  (2.4a)  and  (2.4b),  repeated  here  for  convenience 

xs(tR)  =  x0(tR)  -f  dcos0e(tR)  (B.4) 

ya(tR)  =  y0(tR)  +  d  sin  0e(tR).  (B.5) 

Our  task  now  is  to  solve  Eqs.  (B.1)-(B.5)  for  d,  eliminating  xs(tR )  and  ya{tR).  The 
solution  will  determine  the  two  points  where  the  three  dimensional  ellipsoid  is  intersected 
by  the  line  defined  by  the  intersection  of  the  axial  plane  of  the  effective  beam  and  the 
seabed  plane. 

From  Eq.  (B.3)  we  have 

rT~rR  =  [J/(<:r)  -  **(<r)]2  +  [y/(*r)  -  y*(*R)]2  +  2/(<r) 

-  [*/(**)  -  *i(M]2  -  [y/(*/0  -  y»(*R)]2  -  22(<r) 

=  [z/(*r)  -  *<>(*r)  -  dcos  0e(tR)]2 

+  (y/(*r)  ~  y0(tR.)  -  dsin0e{tR)]2  +  zj(tT) 

-  [*/ (*r)  -  *o(<fl)  -  d cos  0'{tR)]2 

-  [vj^r)  -  Vo{tR)  -  dsin0e{tR)}2  -  z}(tR ) 

=  [z/(*r)  -  z0(<r)]2  +  [y/(<r)  -  !/o(<r)]2  +  2/(<r) 

+  2dcos0e(tR)[xf(tR)  -  xf(tT)] 

-  [*/(**) -  Zo^r)]2  -  [yf(iR.)  -  y0{tR.)]2  -  z2(Ir) 

+  2d  sin  0e(tR)  [yj(tR)  -  yf(tT )] , 

or 

rT  ~  rR  =  2ed  +  /,  (B.6) 

where 

€  =  [*/(0i)“*/(M3 «»*«(**) -[y/(*ii)  +  y/(*r)]  sin  0e(**) 

/  =  lxf(tT)  ~  xo(tR.)]2  +  [yj{h)  -  J/0(<r)]2  +  z){tT) 

-  [*/(<«)  -  xo(<r)]2  -  (j//(<r)  -  yo(tR)]2  -  Z2f{tR). 
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From  Eqs.  (B.l),  (B.2),  (B.4),  and  (B.5)  we  also  have 


rT  +  rR  4  2  rTrR  =  4r2 

=►  2 rTrR  =  4r2  -  r\  -  r2R 

=»  ^rxrft  =  16r4  +  rT  +  rR  +  2rTrR  ~  8r2(rj  +  rR) 

=>  16r4  —  8r2(rj  4  rj-j)  4  rj  4  *"/{  ~  ^rTrR  =  6- 

Therefore, 

16r4  -  8r2(ry  +  rR)  -(-  (r^  -  rR)2  =  0  (B.7) 

From  Eqs.  (B.l),  (B.2),  (B.4),  and  (B.5), 

rT  +  rR  =  [*/(*r)  ~  **(*fl)]2  4  [y/(*r)  -  y»(<*)]2  4  */(<r) 

4  (x/(*/?) _  xa(iR)f  +  [y/(*n)  -  y.(</?)]2  4  ~/(*k) 

=  [xj{tT)  -  x0(tR)  -  dcosde(tR)]2 

4  M*r)  -  y0(iR)  -  d sin  0e(tR))2  4  z)(tT) 

+  (*/(<«)  -  x0(tR)  -  d cos  0e(t r)]2 
4  [y/(<*)  -  y0(*fl)  -  dsin0„(*/*)]2  4  z){tR) 

=  [xsM  -  *o(*k)]2  4  (y/(<r)  -  y0(iR)?  4  z)(tT) 

+  d?  -  2d  cos  0e(tR)  [x/(^k)  +  xj(tT )  -  2x0(tR)} 

4  [*/(**)  ~  Io(^«)]2  4  [y/(<n)  _  yo(^r)]2  +  zj(iR) 

4  d2  -  2d  sin  6e(tn)  [yj(tT)  4  y/(*n)  -  2y0(<«)] , 

or 

rT  +  rfl  =  2d2  —  2 gd  4  2 /i2,  (B.8) 

where 

<7  i  [*/(**)  4  */(*r)  -  2x0(**)]  cos  0e(*R)  4  [y/(**)  4  y/(*r)  -  2y0(*/?)]  sin  0e(tR) 

h2  =  ^  ([*/(*r)  -  ^o(<fl)]2  +  [y/(*r)  -  y0(</i)]2  4  z2(*r) 

+  [*/(</?)  _  ^(^fl)]2  +  [y/(^fl)  -  J/o(</?)]2  +  2r/(^))  • 


Substituting  Eqs.  (B.8)  and  (B.6)  into  (B.7)  we  obtain 

16r4  -  Sr2 (2d2  -  2 gd  +  2 h2)  +  (2 ed  +  f)2  =  0, 
which,  after  rearranging  the  terms,  yields 

4(e2  —  4 r2)d2  +  4(e/  -|-  4 gr2)d  +  16r4  —  16 h2r2  +  f2  =  0. 

The  solution  to  this  quadratic  equation  in  d  yields 

j  —4 (ef  -|-  4 gr2)  ±  ^/l6 (e/  +  4 gr2)2  —  16(e2  —  4r2)(16r4  —  16/i2r2  +  /2) 

8(e2  -  4r2) 

and  after  simplification  we  obtain  the  final  result 

e/  -  gr 2  ±  2r^/\6r4  +  4($2  -  e2  -  4 /i2)r2  +  /2  -f  2e/5  +  4e2/i2 
<*=  2(r2  -  e2)  • 
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